1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Utility functions for the perturbation calculations. 8!> \note 9!> - routines are programmed with spins in mind 10!> but are as of now not tested with them 11!> \par History 12!> 22-08-2002, TCH, started development 13! ************************************************************************************************** 14MODULE qs_p_env_methods 15 USE cp_blacs_env, ONLY: cp_blacs_env_type 16 USE cp_control_types, ONLY: dft_control_type 17 USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,& 18 cp_dbcsr_sm_fm_multiply,& 19 dbcsr_allocate_matrix_set 20 USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,& 21 cp_fm_symm,& 22 cp_fm_triangular_multiply 23 USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose 24 USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,& 25 cp_fm_pool_type,& 26 fm_pool_create_fm,& 27 fm_pool_get_el_struct,& 28 fm_pool_give_back_fm,& 29 fm_pools_create_fm_vect 30 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 31 cp_fm_struct_get,& 32 cp_fm_struct_release,& 33 cp_fm_struct_type 34 USE cp_fm_types, ONLY: cp_fm_create,& 35 cp_fm_get_info,& 36 cp_fm_p_type,& 37 cp_fm_release,& 38 cp_fm_retain,& 39 cp_fm_set_all,& 40 cp_fm_to_fm,& 41 cp_fm_type 42 USE cp_fm_vect, ONLY: cp_fm_vect_copy,& 43 cp_fm_vect_dealloc 44 USE cp_gemm_interface, ONLY: cp_gemm 45 USE cp_log_handling, ONLY: cp_get_default_logger,& 46 cp_logger_type,& 47 cp_to_string 48 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 49 cp_print_key_unit_nr 50 USE cp_para_types, ONLY: cp_para_env_type 51 USE dbcsr_api, ONLY: dbcsr_copy,& 52 dbcsr_p_type,& 53 dbcsr_scale,& 54 dbcsr_set,& 55 dbcsr_type 56 USE hartree_local_methods, ONLY: init_coulomb_local 57 USE hartree_local_types, ONLY: hartree_local_create 58 USE input_constants, ONLY: ot_precond_none 59 USE input_section_types, ONLY: section_vals_get,& 60 section_vals_get_subs_vals,& 61 section_vals_type 62 USE kinds, ONLY: dp 63 USE preconditioner_types, ONLY: init_preconditioner 64 USE qs_energy_types, ONLY: qs_energy_type 65 USE qs_environment_types, ONLY: get_qs_env,& 66 qs_environment_type 67 USE qs_kpp1_env_methods, ONLY: kpp1_calc_k_p_p1,& 68 kpp1_calc_k_p_p1_fdiff,& 69 kpp1_create,& 70 kpp1_did_change 71 USE qs_kpp1_env_types, ONLY: qs_kpp1_env_type 72 USE qs_ks_methods, ONLY: qs_ks_update_qs_env 73 USE qs_ks_types, ONLY: qs_ks_did_change,& 74 qs_ks_env_type 75 USE qs_linres_types, ONLY: linres_control_type 76 USE qs_local_rho_types, ONLY: local_rho_set_create 77 USE qs_matrix_pools, ONLY: mpools_get 78 USE qs_mo_types, ONLY: get_mo_set,& 79 mo_set_p_type 80 USE qs_p_env_types, ONLY: qs_p_env_type 81 USE qs_rho0_methods, ONLY: init_rho0 82 USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals 83 USE qs_rho_methods, ONLY: qs_rho_rebuild,& 84 qs_rho_update_rho 85 USE qs_rho_types, ONLY: qs_rho_create,& 86 qs_rho_get,& 87 qs_rho_type 88 USE string_utilities, ONLY: compress 89#include "./base/base_uses.f90" 90 91 IMPLICIT NONE 92 93 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_p_env_methods' 94 INTEGER, PRIVATE, SAVE :: last_p_env_id = 0 95 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE. 96 97 PRIVATE 98 PUBLIC :: p_env_create, p_env_psi0_changed 99 PUBLIC :: p_op_l1, p_op_l2, p_preortho, p_postortho 100 101CONTAINS 102 103! ************************************************************************************************** 104!> \brief allocates and initializes the perturbation environment (no setup) 105!> \param p_env the environment to initialize 106!> \param qs_env the qs_environment for the system 107!> \param kpp1_env the environment that builds the second order 108!> perturbation kernel 109!> \param p1_option ... 110!> \param psi0d ... 111!> \param orthogonal_orbitals if the orbitals are orthogonal 112!> \param linres_control ... 113!> \par History 114!> 07.2002 created [fawzi] 115!> \author Fawzi Mohamed 116! ************************************************************************************************** 117 SUBROUTINE p_env_create(p_env, qs_env, kpp1_env, p1_option, & 118 psi0d, orthogonal_orbitals, linres_control) 119 120 TYPE(qs_p_env_type), POINTER :: p_env 121 TYPE(qs_environment_type), POINTER :: qs_env 122 TYPE(qs_kpp1_env_type), OPTIONAL, POINTER :: kpp1_env 123 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, & 124 POINTER :: p1_option 125 TYPE(cp_fm_p_type), DIMENSION(:), OPTIONAL, & 126 POINTER :: psi0d 127 LOGICAL, INTENT(in), OPTIONAL :: orthogonal_orbitals 128 TYPE(linres_control_type), OPTIONAL, POINTER :: linres_control 129 130 CHARACTER(len=*), PARAMETER :: routineN = 'p_env_create', routineP = moduleN//':'//routineN 131 132 INTEGER :: handle, n_ao, n_mo, n_spins, natom, spin 133 TYPE(cp_blacs_env_type), POINTER :: blacs_env 134 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools, mo_mo_fm_pools 135 TYPE(cp_fm_type), POINTER :: qs_env_c 136 TYPE(cp_para_env_type), POINTER :: para_env 137 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 138 TYPE(dft_control_type), POINTER :: dft_control 139 140! code 141 142 CALL timeset(routineN, handle) 143 NULLIFY (ao_mo_fm_pools, mo_mo_fm_pools, matrix_s, dft_control, para_env, blacs_env) 144 CALL get_qs_env(qs_env, & 145 matrix_s=matrix_s, & 146 dft_control=dft_control, & 147 para_env=para_env, & 148 blacs_env=blacs_env) 149 150 n_spins = dft_control%nspins 151 152 ALLOCATE (p_env) 153 NULLIFY (p_env%kpp1, & 154 p_env%p1, & 155 p_env%m_epsilon, & 156 p_env%psi0d, & 157 p_env%S_psi0, & 158 p_env%kpp1_env, & 159 p_env%rho1, & 160 p_env%rho1_xc, & 161 p_env%Smo_inv, & 162 p_env%local_rho_set, & 163 p_env%hartree_local, & 164 p_env%PS_psi0, & 165 p_env%preconditioner, & 166 p_env%ev_h0) 167 168 p_env%ref_count = 1 169 last_p_env_id = last_p_env_id + 1 170 p_env%id_nr = last_p_env_id 171 p_env%iter = 0 172 173 p_env%new_preconditioner = .TRUE. 174 p_env%only_energy = .FALSE. 175 p_env%os_valid = .FALSE. 176 p_env%ls_count = 0 177 p_env%delta = 0.0_dp 178 p_env%gnorm = 0.0_dp 179 p_env%gnorm_old = 0.0_dp 180 p_env%etotal = 0.0_dp 181 p_env%gradient = 0.0_dp 182 183 CALL qs_rho_create(p_env%rho1) 184 CALL qs_rho_create(p_env%rho1_xc) 185 186 IF (PRESENT(kpp1_env)) THEN 187 p_env%kpp1_env => kpp1_env 188 ELSE 189 CALL kpp1_create(p_env%kpp1_env) 190 END IF 191 192 IF (PRESENT(p1_option)) THEN 193 p_env%p1 => p1_option 194 ELSE 195 CALL dbcsr_allocate_matrix_set(p_env%p1, n_spins) 196 DO spin = 1, n_spins 197 ALLOCATE (p_env%p1(spin)%matrix) 198 CALL dbcsr_copy(p_env%p1(spin)%matrix, matrix_s(1)%matrix, & 199 name="p_env"//TRIM(ADJUSTL(cp_to_string(p_env%id_nr)))// & 200 "%p1-"//TRIM(ADJUSTL(cp_to_string(spin)))) 201 CALL dbcsr_set(p_env%p1(spin)%matrix, 0.0_dp) 202 END DO 203 END IF 204 205 CALL mpools_get(qs_env%mpools, ao_mo_fm_pools=ao_mo_fm_pools, & 206 mo_mo_fm_pools=mo_mo_fm_pools) 207 208 p_env%n_mo = 0 209 p_env%n_ao = 0 210 DO spin = 1, n_spins 211 IF (PRESENT(psi0d)) THEN 212 CALL cp_fm_get_info(psi0d(spin)%matrix, & 213 ncol_global=n_mo, nrow_global=n_ao) 214 ELSE 215 CALL get_mo_set(qs_env%mos(spin)%mo_set, mo_coeff=qs_env_c) 216 CALL cp_fm_get_info(qs_env_c, & 217 ncol_global=n_mo, nrow_global=n_ao) 218 END IF 219 p_env%n_mo(spin) = n_mo 220 p_env%n_ao(spin) = n_ao 221 END DO 222 223 p_env%orthogonal_orbitals = .FALSE. 224 IF (PRESENT(orthogonal_orbitals)) & 225 p_env%orthogonal_orbitals = orthogonal_orbitals 226 227 CALL fm_pools_create_fm_vect(ao_mo_fm_pools, elements=p_env%S_psi0, & 228 name="p_env"//TRIM(ADJUSTL(cp_to_string(p_env%id_nr)))//"%S_psi0") 229 230 ! alloc m_epsilon 231 CALL fm_pools_create_fm_vect(mo_mo_fm_pools, elements=p_env%m_epsilon, & 232 name="p_env"//TRIM(ADJUSTL(cp_to_string(p_env%id_nr)))// & 233 "%m_epsilon") 234 235 ! alloc Smo_inv 236 IF (.NOT. p_env%orthogonal_orbitals) THEN 237 CALL fm_pools_create_fm_vect(mo_mo_fm_pools, elements=p_env%Smo_inv, & 238 name="p_env"//TRIM(ADJUSTL(cp_to_string(p_env%id_nr)))// & 239 "%Smo_inv") 240 END IF 241 242 IF (PRESENT(psi0d)) THEN 243 IF (ASSOCIATED(psi0d)) THEN 244 CALL cp_fm_vect_copy(psi0d, p_env%psi0d) 245 END IF 246 ELSE IF (.NOT. p_env%orthogonal_orbitals) THEN 247 CALL fm_pools_create_fm_vect(ao_mo_fm_pools, & 248 elements=p_env%psi0d, & 249 name="p_env"//TRIM(ADJUSTL(cp_to_string(p_env%id_nr)))// & 250 "%psi0d") 251 END IF 252 253 !----------------------! 254 ! GAPW initializations ! 255 !----------------------! 256 IF (dft_control%qs_control%gapw) THEN 257 CALL local_rho_set_create(p_env%local_rho_set) 258 CALL allocate_rho_atom_internals(qs_env, p_env%local_rho_set%rho_atom_set) 259 CALL init_rho0(qs_env, dft_control%qs_control%gapw_control, & 260 .TRUE., p_env%local_rho_set) 261 CALL hartree_local_create(p_env%hartree_local) 262 CALL get_qs_env(qs_env=qs_env, natom=natom) 263 CALL init_coulomb_local(p_env%hartree_local, natom) 264 END IF 265 266 !------------------------! 267 ! LINRES initializations ! 268 !------------------------! 269 IF (PRESENT(linres_control)) THEN 270 271 IF (linres_control%preconditioner_type /= ot_precond_none) THEN 272 ! Initialize the preconditioner matrix 273 IF (.NOT. ASSOCIATED(p_env%preconditioner)) THEN 274 275 ALLOCATE (p_env%preconditioner(n_spins)) 276 DO spin = 1, n_spins 277 CALL init_preconditioner(p_env%preconditioner(spin), & 278 para_env=para_env, blacs_env=blacs_env) 279 END DO 280 p_env%os_valid = .FALSE. 281 282 CALL fm_pools_create_fm_vect(ao_mo_fm_pools, elements=p_env%PS_psi0, & 283 name="p_env"//TRIM(ADJUSTL(cp_to_string(p_env%id_nr)))//"%PS_psi0") 284 285 ALLOCATE (p_env%ev_h0(n_spins)) 286 END IF 287 END IF 288 289 END IF 290 291 CALL timestop(handle) 292 293 END SUBROUTINE p_env_create 294 295! ************************************************************************************************** 296!> \brief checks that the intenal storage is allocated, and allocs it if needed 297!> \param p_env the environment to check 298!> \param qs_env the qs environment this p_env lives in 299!> \par History 300!> 12.2002 created [fawzi] 301!> \author Fawzi Mohamed 302!> \note 303!> private routine 304! ************************************************************************************************** 305 SUBROUTINE p_env_check_i_alloc(p_env, qs_env) 306 TYPE(qs_p_env_type), POINTER :: p_env 307 TYPE(qs_environment_type), POINTER :: qs_env 308 309 CHARACTER(len=*), PARAMETER :: routineN = 'p_env_check_i_alloc', & 310 routineP = moduleN//':'//routineN 311 312 CHARACTER(len=25) :: name 313 INTEGER :: handle, ispin, nspins 314 LOGICAL :: gapw_xc 315 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 316 TYPE(dft_control_type), POINTER :: dft_control 317 318 CALL timeset(routineN, handle) 319 320 NULLIFY (dft_control, matrix_s) 321 322 CPASSERT(ASSOCIATED(p_env)) 323 CPASSERT(p_env%ref_count > 0) 324 CALL get_qs_env(qs_env, dft_control=dft_control) 325 gapw_xc = dft_control%qs_control%gapw_xc 326 IF (.NOT. ASSOCIATED(p_env%kpp1)) THEN 327 CALL get_qs_env(qs_env, matrix_s=matrix_s) 328 nspins = dft_control%nspins 329 330 CALL dbcsr_allocate_matrix_set(p_env%kpp1, nspins) 331 name = "p_env"//cp_to_string(p_env%id_nr)//"%kpp1-" 332 CALL compress(name, full=.TRUE.) 333 DO ispin = 1, nspins 334 ALLOCATE (p_env%kpp1(ispin)%matrix) 335 CALL dbcsr_copy(p_env%kpp1(ispin)%matrix, matrix_s(1)%matrix, & 336 name=TRIM(name)//ADJUSTL(cp_to_string(ispin))) 337 END DO 338 339 CALL qs_rho_rebuild(p_env%rho1, qs_env=qs_env) 340 IF (gapw_xc) THEN 341 CALL qs_rho_rebuild(p_env%rho1_xc, qs_env=qs_env) 342 END IF 343 344 END IF 345 CALL timestop(handle) 346 END SUBROUTINE p_env_check_i_alloc 347 348! ************************************************************************************************** 349!> \brief To be called after the value of psi0 has changed. 350!> Recalculates the quantities S_psi0 and m_epsilon. 351!> \param p_env the perturbation environment to set 352!> \param qs_env ... 353!> \param psi0 the value of psi0, if not given defaults to the qs_env mos 354!> \param Hrho_psi0d is given, then the partial result Hrho_psi0d is stored in 355!> that vector 356!> \par History 357!> 07.2002 created [fawzi] 358!> \author Fawzi Mohamed 359! ************************************************************************************************** 360 SUBROUTINE p_env_psi0_changed(p_env, qs_env, psi0, Hrho_psi0d) 361 362 TYPE(qs_p_env_type), POINTER :: p_env 363 TYPE(qs_environment_type), POINTER :: qs_env 364 TYPE(cp_fm_p_type), DIMENSION(:), OPTIONAL, & 365 POINTER :: psi0 366 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(inout), & 367 OPTIONAL :: Hrho_psi0d 368 369 CHARACTER(len=*), PARAMETER :: routineN = 'p_env_psi0_changed', & 370 routineP = moduleN//':'//routineN 371 372 INTEGER :: handle, iounit, lfomo, n_spins, nmo, spin 373 LOGICAL :: was_present 374 REAL(KIND=dp) :: maxocc 375 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: my_psi0 376 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_mo_fm_pools 377 TYPE(cp_logger_type), POINTER :: logger 378 TYPE(cp_para_env_type), POINTER :: para_env 379 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao 380 TYPE(dft_control_type), POINTER :: dft_control 381 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 382 TYPE(qs_energy_type), POINTER :: energy 383 TYPE(qs_ks_env_type), POINTER :: ks_env 384 TYPE(qs_rho_type), POINTER :: rho 385 TYPE(section_vals_type), POINTER :: input, lr_section 386 387 CALL timeset(routineN, handle) 388 389 NULLIFY (ao_mo_fm_pools, mos, my_psi0, matrix_s, mos, para_env, ks_env, rho, & 390 logger, input, lr_section, energy, matrix_ks, dft_control, rho_ao) 391 logger => cp_get_default_logger() 392 393 CPASSERT(ASSOCIATED(p_env)) 394 CPASSERT(p_env%ref_count > 0) 395 396 CALL get_qs_env(qs_env, & 397 ks_env=ks_env, & 398 mos=mos, & 399 matrix_s=matrix_s, & 400 matrix_ks=matrix_ks, & 401 para_env=para_env, & 402 rho=rho, & 403 input=input, & 404 energy=energy, & 405 dft_control=dft_control) 406 407 CALL qs_rho_get(rho, rho_ao=rho_ao) 408 409 n_spins = dft_control%nspins 410 p_env%iter = p_env%iter + 1 411 CALL mpools_get(qs_env%mpools, & 412 ao_mo_fm_pools=ao_mo_fm_pools) 413 ! def my_psi0 414 IF (PRESENT(psi0)) THEN 415 CALL cp_fm_vect_copy(psi0, my_psi0) 416 ELSE 417 ALLOCATE (my_psi0(n_spins)) 418 DO spin = 1, n_spins 419 NULLIFY (my_psi0(spin)%matrix) 420 CALL get_mo_set(mos(spin)%mo_set, & 421 mo_coeff=my_psi0(spin)%matrix) 422 CALL cp_fm_retain(my_psi0(spin)%matrix) 423 END DO 424 END IF 425 426 lr_section => section_vals_get_subs_vals(input, "PROPERTIES%LINRES") 427 ! def psi0d 428 IF (p_env%orthogonal_orbitals) THEN 429 IF (ASSOCIATED(p_env%psi0d)) THEN 430 CALL cp_fm_vect_dealloc(p_env%psi0d) 431 END IF 432 p_env%psi0d => my_psi0 433 ELSE 434 435 DO spin = 1, n_spins 436 ! m_epsilon=cholesky_decomposition(my_psi0^T S my_psi0)^-1 437 ! could be optimized by combining next two calls 438 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, & 439 my_psi0(spin)%matrix, & 440 p_env%S_psi0(spin)%matrix, & 441 ncol=p_env%n_mo(spin), alpha=1.0_dp) 442 CALL cp_gemm(transa='T', transb='N', n=p_env%n_mo(spin), & 443 m=p_env%n_mo(spin), k=p_env%n_ao(spin), alpha=1.0_dp, & 444 matrix_a=my_psi0(spin)%matrix, & 445 matrix_b=p_env%S_psi0(spin)%matrix, & 446 beta=0.0_dp, matrix_c=p_env%m_epsilon(spin)%matrix) 447 CALL cp_fm_cholesky_decompose(p_env%m_epsilon(spin)%matrix, & 448 n=p_env%n_mo(spin)) 449 450 ! Smo_inv= (my_psi0^T S my_psi0)^-1 451 CALL cp_fm_set_all(p_env%Smo_inv(spin)%matrix, 0.0_dp, 1.0_dp) 452 ! faster using cp_fm_cholesky_invert ? 453 CALL cp_fm_triangular_multiply( & 454 triangular_matrix=p_env%m_epsilon(spin)%matrix, & 455 matrix_b=p_env%Smo_inv(spin)%matrix, side='R', & 456 invert_tr=.TRUE., n_rows=p_env%n_mo(spin), & 457 n_cols=p_env%n_mo(spin)) 458 CALL cp_fm_triangular_multiply( & 459 triangular_matrix=p_env%m_epsilon(spin)%matrix, & 460 matrix_b=p_env%Smo_inv(spin)%matrix, side='R', & 461 transpose_tr=.TRUE., & 462 invert_tr=.TRUE., n_rows=p_env%n_mo(spin), & 463 n_cols=p_env%n_mo(spin)) 464 465 ! psi0d=my_psi0 (my_psi0^T S my_psi0)^-1 466 ! faster using cp_fm_cholesky_invert ? 467 CALL cp_fm_to_fm(my_psi0(spin)%matrix, & 468 p_env%psi0d(spin)%matrix) 469 CALL cp_fm_triangular_multiply( & 470 triangular_matrix=p_env%m_epsilon(spin)%matrix, & 471 matrix_b=p_env%psi0d(spin)%matrix, side='R', & 472 invert_tr=.TRUE., n_rows=p_env%n_ao(spin), & 473 n_cols=p_env%n_mo(spin)) 474 CALL cp_fm_triangular_multiply( & 475 triangular_matrix=p_env%m_epsilon(spin)%matrix, & 476 matrix_b=p_env%psi0d(spin)%matrix, side='R', & 477 transpose_tr=.TRUE., & 478 invert_tr=.TRUE., n_rows=p_env%n_ao(spin), & 479 n_cols=p_env%n_mo(spin)) 480 481 ! updates P 482 CALL get_mo_set(mos(spin)%mo_set, lfomo=lfomo, & 483 nmo=nmo, maxocc=maxocc) 484 IF (lfomo > nmo) THEN 485 CALL dbcsr_set(rho_ao(spin)%matrix, 0.0_dp) 486 CALL cp_dbcsr_plus_fm_fm_t(rho_ao(spin)%matrix, & 487 matrix_v=my_psi0(spin)%matrix, & 488 matrix_g=p_env%psi0d(spin)%matrix, & 489 ncol=p_env%n_mo(spin)) 490 CALL dbcsr_scale(rho_ao(spin)%matrix, alpha_scalar=maxocc) 491 ELSE 492 CPABORT("symmetrized onesided smearing to do") 493 END IF 494 END DO 495 496 ! updates rho 497 CALL qs_rho_update_rho(rho_struct=rho, qs_env=qs_env) 498 499 ! tells ks_env that p changed 500 CALL qs_ks_did_change(ks_env=ks_env, rho_changed=.TRUE.) 501 502 END IF 503 504 ! updates K (if necessary) 505 CALL qs_ks_update_qs_env(qs_env) 506 iounit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", & 507 extension=".linresLog") 508 IF (iounit > 0) THEN 509 CALL section_vals_get(lr_section, explicit=was_present) 510 IF (was_present) THEN 511 WRITE (UNIT=iounit, FMT="(/,(T3,A,T55,F25.14))") & 512 "Total energy ground state: ", energy%total 513 END IF 514 END IF 515 CALL cp_print_key_finished_output(iounit, logger, lr_section, & 516 "PRINT%PROGRAM_RUN_INFO") 517 !-----------------------------------------------------------------------| 518 ! calculates | 519 ! m_epsilon = - psi0d^T times K times psi0d | 520 ! = - [K times psi0d]^T times psi0d (because K is symmetric) | 521 !-----------------------------------------------------------------------| 522 DO spin = 1, n_spins 523 ! S_psi0 = k times psi0d 524 CALL cp_dbcsr_sm_fm_multiply(matrix_ks(spin)%matrix, & 525 p_env%psi0d(spin)%matrix, & 526 p_env%S_psi0(spin)%matrix, p_env%n_mo(spin)) 527 IF (PRESENT(Hrho_psi0d)) THEN 528 CALL cp_fm_scale_and_add(alpha=0.0_dp, matrix_a=Hrho_psi0d(spin)%matrix, & 529 beta=1.0_dp, matrix_b=p_env%S_psi0(spin)%matrix) 530 END IF 531 ! m_epsilon = -1 times S_psi0^T times psi0d 532 CALL cp_gemm('T', 'N', & 533 p_env%n_mo(spin), p_env%n_mo(spin), p_env%n_ao(spin), & 534 -1.0_dp, p_env%S_psi0(spin)%matrix, p_env%psi0d(spin)%matrix, & 535 0.0_dp, p_env%m_epsilon(spin)%matrix) 536 END DO 537 538 !----------------------------------| 539 ! calculates S_psi0 = S * my_psi0 | 540 !----------------------------------| 541 ! calculating this reduces the mat mult without storing a full aoxao 542 ! matrix (for P). If nspin>1 you might consider calculating it on the 543 ! fly to spare some memory 544 CALL get_qs_env(qs_env, matrix_s=matrix_s) 545 DO spin = 1, n_spins 546 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, & 547 my_psi0(spin)%matrix, & 548 p_env%S_psi0(spin)%matrix, & 549 p_env%n_mo(spin)) 550 END DO 551 552 ! releases my_psi0 553 IF (p_env%orthogonal_orbitals) THEN 554 NULLIFY (my_psi0) 555 ELSE 556 CALL cp_fm_vect_dealloc(my_psi0) 557 END IF 558 559 ! tells kpp1_env about the change of psi0 560 CALL kpp1_did_change(p_env%kpp1_env, psi0_changed=.TRUE.) 561 562 CALL timestop(handle) 563 564 END SUBROUTINE p_env_psi0_changed 565 566! ************************************************************************************************** 567!> \brief Evaluates Fv (S_mo)^-1 - Sv(epsilon) and stores it in res 568!> \param p_env perturbation calculation environment 569!> \param qs_env the qs_env that is perturbed by this p_env 570!> \param v the matrix to operate on 571!> \param res the result 572!> \par History 573!> 10.2002, TCH, extracted single spin calculation 574!> \author Thomas Chassaing 575! ************************************************************************************************** 576 SUBROUTINE p_op_l1(p_env, qs_env, v, res) 577 578 ! argument 579 TYPE(qs_p_env_type), POINTER :: p_env 580 TYPE(qs_environment_type), POINTER :: qs_env 581 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(in) :: v 582 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(inout) :: res 583 584 CHARACTER(len=*), PARAMETER :: routineN = 'p_op_l1', routineP = moduleN//':'//routineN 585 586 INTEGER :: n_spins, spin 587 TYPE(dft_control_type), POINTER :: dft_control 588 589 CPASSERT(ASSOCIATED(p_env)) 590 CPASSERT(p_env%ref_count > 0) 591 NULLIFY (dft_control) 592 CALL get_qs_env(qs_env, dft_control=dft_control) 593 594 n_spins = dft_control%nspins 595 DO spin = 1, n_spins 596 CALL p_op_l1_spin(p_env, qs_env, spin, v(spin)%matrix, & 597 res(spin)%matrix) 598 END DO 599 600 END SUBROUTINE p_op_l1 601 602! ************************************************************************************************** 603!> \brief Evaluates Fv (S_mo)^-1 - Sv(epsilon) and stores it in res 604!> for a given spin 605!> \param p_env perturbation calculation environment 606!> \param qs_env the qs_env that is perturbed by this p_env 607!> \param spin the spin to calculate (1 or 2 normally) 608!> \param v the matrix to operate on 609!> \param res the result 610!> \par History 611!> 10.2002, TCH, created 612!> \author Thomas Chassaing 613!> \note 614!> Same as p_op_l1 but takes a spin as additional argument. 615! ************************************************************************************************** 616 SUBROUTINE p_op_l1_spin(p_env, qs_env, spin, v, res) 617 618 ! argument 619 TYPE(qs_p_env_type), POINTER :: p_env 620 TYPE(qs_environment_type), POINTER :: qs_env 621 INTEGER, INTENT(IN) :: spin 622 TYPE(cp_fm_type), POINTER :: v, res 623 624 CHARACTER(len=*), PARAMETER :: routineN = 'p_op_l1_spin', routineP = moduleN//':'//routineN 625 626 INTEGER :: handle, ncol 627 TYPE(cp_fm_type), POINTER :: tmp 628 TYPE(cp_para_env_type), POINTER :: para_env 629 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s 630 TYPE(dbcsr_type), POINTER :: k_p 631 TYPE(dft_control_type), POINTER :: dft_control 632 633 CALL timeset(routineN, handle) 634 NULLIFY (tmp, matrix_ks, matrix_s, dft_control) 635 636 CALL get_qs_env(qs_env, & 637 para_env=para_env, & 638 matrix_s=matrix_s, & 639 matrix_ks=matrix_ks, & 640 dft_control=dft_control) 641 642 CPASSERT(ASSOCIATED(p_env)) 643 CPASSERT(p_env%ref_count > 0) 644 CPASSERT(0 < spin) 645 CPASSERT(spin <= dft_control%nspins) 646 647 CALL cp_fm_create(tmp, v%matrix_struct) 648 649 k_p => matrix_ks(spin)%matrix 650 CALL cp_fm_get_info(v, ncol_global=ncol) 651 652 IF (p_env%orthogonal_orbitals) THEN 653 CALL cp_dbcsr_sm_fm_multiply(k_p, v, res, ncol) 654 ELSE 655 CALL cp_dbcsr_sm_fm_multiply(k_p, v, tmp, ncol) 656 CALL cp_fm_symm('R', 'U', p_env%n_ao(spin), p_env%n_mo(spin), 1.0_dp, & 657 p_env%Smo_inv(spin)%matrix, tmp, 0.0_dp, res) 658 END IF 659 660 CALL cp_fm_symm('R', 'U', p_env%n_ao(spin), p_env%n_mo(spin), 1.0_dp, & 661 p_env%m_epsilon(spin)%matrix, v, 0.0_dp, tmp) 662 CALL cp_dbcsr_sm_fm_multiply(matrix_s(1)%matrix, tmp, & 663 res, p_env%n_mo(spin), alpha=1.0_dp, beta=1.0_dp) 664 665 CALL cp_fm_release(tmp) 666 CALL timestop(handle) 667 668 END SUBROUTINE p_op_l1_spin 669 670! ************************************************************************************************** 671!> \brief evaluates res = alpha kpp1(v)*psi0 + beta res 672!> with kpp1 evaluated with p=qs_env%rho%rho_ao, p1=p1 673!> \param p_env the perturbation environment 674!> \param qs_env the qs_env that is perturbed by this p_env 675!> \param p1 direction in which evaluate the second derivative 676!> \param res place where to store the result 677!> \param alpha scale factor of the result (defaults to 1.0) 678!> \param beta scale factor of the old values (defaults to 0.0) 679!> \par History 680!> 02.09.2002 adapted for new qs_p_env_type (TC) 681!> 03.2003 extended for p1 not taken from v (TC) 682!> \author fawzi 683!> \note 684!> qs_env%rho must be up to date 685!> it would be better to pass rho1, not p1 686! ************************************************************************************************** 687 SUBROUTINE p_op_l2(p_env, qs_env, p1, res, alpha, beta) 688 689 TYPE(qs_p_env_type), POINTER :: p_env 690 TYPE(qs_environment_type), POINTER :: qs_env 691 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: p1 692 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(INOUT) :: res 693 REAL(KIND=dp), INTENT(in), OPTIONAL :: alpha, beta 694 695 CHARACTER(len=*), PARAMETER :: routineN = 'p_op_l2', routineP = moduleN//':'//routineN 696 LOGICAL, PARAMETER :: fdiff = .FALSE. 697 698 INTEGER :: handle, ispin, n_spins 699 INTEGER, SAVE :: iter = 0 700 LOGICAL :: gapw, gapw_xc 701 REAL(KIND=dp) :: my_alpha, my_beta 702 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho1_ao 703 TYPE(dft_control_type), POINTER :: dft_control 704 TYPE(qs_rho_type), POINTER :: rho 705 706 CALL timeset(routineN, handle) 707 NULLIFY (dft_control, rho, rho1_ao) 708 709 CALL get_qs_env(qs_env, dft_control=dft_control, rho=rho) 710 711 gapw = dft_control%qs_control%gapw 712 gapw_xc = dft_control%qs_control%gapw_xc 713 my_alpha = 1.0_dp 714 IF (PRESENT(alpha)) my_alpha = alpha 715 my_beta = 0.0_dp 716 IF (PRESENT(beta)) my_beta = beta 717 718 iter = iter + 1 719 720 CPASSERT(ASSOCIATED(p_env)) 721 CPASSERT(p_env%ref_count > 0) 722 723 CALL p_env_check_i_alloc(p_env, qs_env=qs_env) 724 n_spins = dft_control%nspins 725 726 CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao) 727 DO ispin = 1, SIZE(p1) 728 ! hack to avoid crashes in ep regs 729 IF (.NOT. ASSOCIATED(rho1_ao(ispin)%matrix, p1(ispin)%matrix)) THEN 730 CALL dbcsr_copy(rho1_ao(ispin)%matrix, p1(ispin)%matrix) 731 ENDIF 732 ENDDO 733 734 IF (ASSOCIATED(p_env%local_rho_set)) THEN 735 CALL qs_rho_update_rho(rho_struct=p_env%rho1, qs_env=qs_env, local_rho_set=p_env%local_rho_set) 736 ELSE 737 CALL qs_rho_update_rho(rho_struct=p_env%rho1, qs_env=qs_env) 738 END IF 739 740 IF (fdiff) THEN 741 CALL kpp1_calc_k_p_p1_fdiff(qs_env=qs_env, & 742 k_p_p1=p_env%kpp1, rho=rho, rho1=p_env%rho1) 743 ELSE 744 CALL kpp1_calc_k_p_p1(kpp1_env=p_env%kpp1_env, p_env=p_env, qs_env=qs_env, & 745 k_p_p1=p_env%kpp1, rho=rho, rho1=p_env%rho1, rho1_xc=p_env%rho1) 746 END IF 747 748 DO ispin = 1, n_spins 749 CALL cp_dbcsr_sm_fm_multiply(p_env%kpp1(ispin)%matrix, & 750 p_env%psi0d(ispin)%matrix, res(ispin)%matrix, & 751 ncol=p_env%n_mo(ispin), & 752 alpha=my_alpha, beta=my_beta) 753 END DO 754 755 CALL timestop(handle) 756 757 END SUBROUTINE p_op_l2 758 759! ************************************************************************************************** 760!> \brief does a preorthogonalization of the given matrix: 761!> v = (I-PS)v 762!> \param p_env the perturbation environment 763!> \param qs_env the qs_env that is perturbed by this p_env 764!> \param v matrix to orthogonalize 765!> \param n_cols the number of columns of C to multiply (defaults to size(v,2)) 766!> \par History 767!> 02.09.2002 adapted for new qs_p_env_type (TC) 768!> \author Fawzi Mohamed 769! ************************************************************************************************** 770 SUBROUTINE p_preortho(p_env, qs_env, v, n_cols) 771 772 TYPE(qs_p_env_type), POINTER :: p_env 773 TYPE(qs_environment_type), POINTER :: qs_env 774 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(inout) :: v 775 INTEGER, DIMENSION(:), INTENT(in), OPTIONAL :: n_cols 776 777 CHARACTER(len=*), PARAMETER :: routineN = 'p_preortho', routineP = moduleN//':'//routineN 778 779 INTEGER :: cols, handle, max_cols, maxnmo, n_spins, & 780 nmo2, spin, v_cols, v_rows 781 TYPE(cp_fm_pool_type), POINTER :: maxmo_maxmo_fm_pool 782 TYPE(cp_fm_struct_type), POINTER :: maxmo_maxmo_fmstruct, tmp_fmstruct 783 TYPE(cp_fm_type), POINTER :: tmp_matrix 784 TYPE(dft_control_type), POINTER :: dft_control 785 786 CALL timeset(routineN, handle) 787 788 NULLIFY (tmp_matrix, maxmo_maxmo_fm_pool, maxmo_maxmo_fmstruct, tmp_fmstruct, & 789 dft_control) 790 791 CPASSERT(ASSOCIATED(p_env)) 792 CPASSERT(p_env%ref_count > 0) 793 794 CALL get_qs_env(qs_env, dft_control=dft_control) 795 CALL mpools_get(qs_env%mpools, maxmo_maxmo_fm_pool=maxmo_maxmo_fm_pool) 796 n_spins = dft_control%nspins 797 maxmo_maxmo_fmstruct => fm_pool_get_el_struct(maxmo_maxmo_fm_pool) 798 CALL cp_fm_struct_get(maxmo_maxmo_fmstruct, nrow_global=nmo2, ncol_global=maxnmo) 799 CPASSERT(SIZE(v) >= n_spins) 800 ! alloc tmp storage 801 IF (PRESENT(n_cols)) THEN 802 max_cols = MAXVAL(n_cols(1:n_spins)) 803 ELSE 804 max_cols = 0 805 DO spin = 1, n_spins 806 CALL cp_fm_get_info(v(spin)%matrix, ncol_global=v_cols) 807 max_cols = MAX(max_cols, v_cols) 808 END DO 809 END IF 810 IF (max_cols <= nmo2) THEN 811 CALL fm_pool_create_fm(maxmo_maxmo_fm_pool, tmp_matrix) 812 ELSE 813 CALL cp_fm_struct_create(tmp_fmstruct, nrow_global=max_cols, & 814 ncol_global=maxnmo, template_fmstruct=maxmo_maxmo_fmstruct) 815 CALL cp_fm_create(tmp_matrix, matrix_struct=tmp_fmstruct) 816 CALL cp_fm_struct_release(tmp_fmstruct) 817 END IF 818 819 DO spin = 1, n_spins 820 821 CALL cp_fm_get_info(v(spin)%matrix, & 822 nrow_global=v_rows, ncol_global=v_cols) 823 CPASSERT(v_rows >= p_env%n_ao(spin)) 824 cols = v_cols 825 IF (PRESENT(n_cols)) THEN 826 CPASSERT(n_cols(spin) <= cols) 827 cols = n_cols(spin) 828 END IF 829 CPASSERT(cols <= max_cols) 830 831 ! tmp_matrix = v^T (S psi0) 832 CALL cp_gemm(transa='T', transb='N', m=cols, n=p_env%n_mo(spin), & 833 k=p_env%n_ao(spin), alpha=1.0_dp, matrix_a=v(spin)%matrix, & 834 matrix_b=p_env%S_psi0(spin)%matrix, beta=0.0_dp, & 835 matrix_c=tmp_matrix) 836 ! v = v - psi0d tmp_matrix^T = v - psi0d psi0^T S v 837 CALL cp_gemm(transa='N', transb='T', m=p_env%n_ao(spin), n=cols, & 838 k=p_env%n_mo(spin), alpha=-1.0_dp, & 839 matrix_a=p_env%psi0d(spin)%matrix, matrix_b=tmp_matrix, & 840 beta=1.0_dp, matrix_c=v(spin)%matrix) 841 842 END DO 843 844 IF (max_cols <= nmo2) THEN 845 CALL fm_pool_give_back_fm(maxmo_maxmo_fm_pool, tmp_matrix) 846 ELSE 847 CALL cp_fm_release(tmp_matrix) 848 END IF 849 850 CALL timestop(handle) 851 852 END SUBROUTINE p_preortho 853 854! ************************************************************************************************** 855!> \brief does a postorthogonalization on the given matrix vector: 856!> v = (I-SP) v 857!> \param p_env the perturbation environment 858!> \param qs_env the qs_env that is perturbed by this p_env 859!> \param v matrix to orthogonalize 860!> \param n_cols the number of columns of C to multiply (defaults to size(v,2)) 861!> \par History 862!> 07.2002 created [fawzi] 863!> \author Fawzi Mohamed 864! ************************************************************************************************** 865 SUBROUTINE p_postortho(p_env, qs_env, v, n_cols) 866 867 TYPE(qs_p_env_type), POINTER :: p_env 868 TYPE(qs_environment_type), POINTER :: qs_env 869 TYPE(cp_fm_p_type), DIMENSION(:), INTENT(inout) :: v 870 INTEGER, DIMENSION(:), INTENT(in), OPTIONAL :: n_cols 871 872 CHARACTER(len=*), PARAMETER :: routineN = 'p_postortho', routineP = moduleN//':'//routineN 873 874 INTEGER :: cols, handle, max_cols, maxnmo, n_spins, & 875 nmo2, spin, v_cols, v_rows 876 TYPE(cp_fm_pool_type), POINTER :: maxmo_maxmo_fm_pool 877 TYPE(cp_fm_struct_type), POINTER :: maxmo_maxmo_fmstruct, tmp_fmstruct 878 TYPE(cp_fm_type), POINTER :: tmp_matrix 879 TYPE(dft_control_type), POINTER :: dft_control 880 881 CALL timeset(routineN, handle) 882 883 NULLIFY (tmp_matrix, maxmo_maxmo_fm_pool, maxmo_maxmo_fmstruct, tmp_fmstruct, & 884 dft_control) 885 886 CPASSERT(ASSOCIATED(p_env)) 887 CPASSERT(p_env%ref_count > 0) 888 889 CALL get_qs_env(qs_env, dft_control=dft_control) 890 CALL mpools_get(qs_env%mpools, maxmo_maxmo_fm_pool=maxmo_maxmo_fm_pool) 891 n_spins = dft_control%nspins 892 maxmo_maxmo_fmstruct => fm_pool_get_el_struct(maxmo_maxmo_fm_pool) 893 CALL cp_fm_struct_get(maxmo_maxmo_fmstruct, nrow_global=nmo2, ncol_global=maxnmo) 894 CPASSERT(SIZE(v) >= n_spins) 895 ! alloc tmp storage 896 IF (PRESENT(n_cols)) THEN 897 max_cols = MAXVAL(n_cols(1:n_spins)) 898 ELSE 899 max_cols = 0 900 DO spin = 1, n_spins 901 CALL cp_fm_get_info(v(spin)%matrix, ncol_global=v_cols) 902 max_cols = MAX(max_cols, v_cols) 903 END DO 904 END IF 905 IF (max_cols <= nmo2) THEN 906 CALL fm_pool_create_fm(maxmo_maxmo_fm_pool, tmp_matrix) 907 ELSE 908 CALL cp_fm_struct_create(tmp_fmstruct, nrow_global=max_cols, & 909 ncol_global=maxnmo, template_fmstruct=maxmo_maxmo_fmstruct) 910 CALL cp_fm_create(tmp_matrix, matrix_struct=tmp_fmstruct) 911 CALL cp_fm_struct_release(tmp_fmstruct) 912 END IF 913 914 DO spin = 1, n_spins 915 916 CALL cp_fm_get_info(v(spin)%matrix, & 917 nrow_global=v_rows, ncol_global=v_cols) 918 CPASSERT(v_rows >= p_env%n_ao(spin)) 919 cols = v_cols 920 IF (PRESENT(n_cols)) THEN 921 CPASSERT(n_cols(spin) <= cols) 922 cols = n_cols(spin) 923 END IF 924 CPASSERT(cols <= max_cols) 925 926 ! tmp_matrix = v^T psi0d 927 CALL cp_gemm(transa='T', transb='N', m=cols, n=p_env%n_mo(spin), & 928 k=p_env%n_ao(spin), alpha=1.0_dp, matrix_a=v(spin)%matrix, & 929 matrix_b=p_env%psi0d(spin)%matrix, beta=0.0_dp, & 930 matrix_c=tmp_matrix) 931 ! v = v - (S psi0) tmp_matrix^T = v - S psi0 psi0d^T v 932 CALL cp_gemm(transa='N', transb='T', m=p_env%n_ao(spin), n=cols, & 933 k=p_env%n_mo(spin), alpha=-1.0_dp, & 934 matrix_a=p_env%S_psi0(spin)%matrix, matrix_b=tmp_matrix, & 935 beta=1.0_dp, matrix_c=v(spin)%matrix) 936 937 END DO 938 939 IF (max_cols <= nmo2) THEN 940 CALL fm_pool_give_back_fm(maxmo_maxmo_fm_pool, tmp_matrix) 941 ELSE 942 CALL cp_fm_release(tmp_matrix) 943 END IF 944 945 CALL timestop(handle) 946 947 END SUBROUTINE p_postortho 948 949END MODULE qs_p_env_methods 950