1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Utility subroutine for qs energy calculation 8!> \par History 9!> none 10!> \author MK (29.10.2002) 11! ************************************************************************************************** 12MODULE qs_energy_utils 13 USE atprop_types, ONLY: atprop_array_add,& 14 atprop_type 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_control_utils, ONLY: read_ddapc_section 17 USE cp_external_control, ONLY: external_control 18 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 19 cp_fm_struct_release,& 20 cp_fm_struct_type 21 USE cp_fm_types, ONLY: cp_fm_create,& 22 cp_fm_p_type,& 23 cp_fm_release,& 24 cp_fm_type 25 USE dbcsr_api, ONLY: & 26 dbcsr_add, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, & 27 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & 28 dbcsr_p_type, dbcsr_set, dbcsr_type 29 USE et_coupling, ONLY: calc_et_coupling 30 USE et_coupling_proj, ONLY: calc_et_coupling_proj 31 USE input_section_types, ONLY: section_vals_get,& 32 section_vals_get_subs_vals,& 33 section_vals_type 34 USE kinds, ONLY: dp 35 USE kpoint_methods, ONLY: kpoint_density_matrices,& 36 kpoint_density_transform 37 USE kpoint_types, ONLY: kpoint_type 38 USE mp2, ONLY: mp2_main 39 USE pw_env_types, ONLY: pw_env_type 40 USE pw_types, ONLY: pw_p_type 41 USE qs_energy_types, ONLY: qs_energy_type 42 USE qs_environment_types, ONLY: get_qs_env,& 43 qs_environment_type 44 USE qs_integrate_potential, ONLY: integrate_v_core_rspace 45 USE qs_ks_methods, ONLY: calculate_w_matrix,& 46 calculate_w_matrix_ot,& 47 qs_ks_update_qs_env 48 USE qs_linres_module, ONLY: linres_calculation_low 49 USE qs_mo_types, ONLY: get_mo_set,& 50 mo_set_p_type,& 51 mo_set_type 52 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 53 USE qs_rho_types, ONLY: qs_rho_get,& 54 qs_rho_type 55 USE qs_scf, ONLY: scf 56 USE qs_tddfpt2_methods, ONLY: tddfpt 57 USE scf_control_types, ONLY: scf_control_type 58 USE xas_methods, ONLY: xas 59 USE xas_tdp_methods, ONLY: xas_tdp 60#include "./base/base_uses.f90" 61 62 IMPLICIT NONE 63 64 PRIVATE 65 66! *** Global parameters *** 67 68 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_energy_utils' 69 70 PUBLIC :: qs_energies_compute_matrix_w, & 71 qs_energies_properties, & 72 qs_energies_mp2 73 74CONTAINS 75 76! ************************************************************************************************** 77!> \brief Refactoring of qs_energies_scf. Moves computation of matrix_w 78!> into separate subroutine 79!> \param qs_env ... 80!> \param calc_forces ... 81!> \par History 82!> 05.2013 created [Florian Schiffmann] 83! ************************************************************************************************** 84 85 SUBROUTINE qs_energies_compute_matrix_w(qs_env, calc_forces) 86 TYPE(qs_environment_type), POINTER :: qs_env 87 LOGICAL, INTENT(IN) :: calc_forces 88 89 CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_compute_matrix_w', & 90 routineP = moduleN//':'//routineN 91 92 INTEGER :: handle, is, ispin, nao, nspin 93 LOGICAL :: do_kpoints, has_unit_metric 94 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: fmwork 95 TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct 96 TYPE(cp_fm_type), POINTER :: mo_coeff 97 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, matrix_w, & 98 matrix_w_mp2, mo_derivs, rho_ao 99 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, matrix_w_kp 100 TYPE(dft_control_type), POINTER :: dft_control 101 TYPE(kpoint_type), POINTER :: kpoints 102 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 103 TYPE(mo_set_type), POINTER :: mo_set 104 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 105 POINTER :: sab_nl 106 TYPE(qs_rho_type), POINTER :: rho 107 TYPE(scf_control_type), POINTER :: scf_control 108 109 CALL timeset(routineN, handle) 110 111 ! if calculate forces, time to compute the w matrix 112 CALL get_qs_env(qs_env, has_unit_metric=has_unit_metric) 113 114 IF (calc_forces .AND. .NOT. has_unit_metric) THEN 115 CALL get_qs_env(qs_env, do_kpoints=do_kpoints) 116 117 IF (do_kpoints) THEN 118 119 CALL get_qs_env(qs_env, & 120 matrix_w_kp=matrix_w_kp, & 121 matrix_s_kp=matrix_s_kp, & 122 sab_orb=sab_nl, & 123 mos=mos, & 124 kpoints=kpoints) 125 126 CALL get_mo_set(mos(1)%mo_set, mo_coeff=mo_coeff, nao=nao) 127 CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, nrow_global=nao, ncol_global=nao, & 128 template_fmstruct=mo_coeff%matrix_struct) 129 130 ALLOCATE (fmwork(2)) 131 DO is = 1, SIZE(fmwork) 132 NULLIFY (fmwork(is)%matrix) 133 CALL cp_fm_create(fmwork(is)%matrix, matrix_struct=ao_ao_fmstruct) 134 END DO 135 CALL cp_fm_struct_release(ao_ao_fmstruct) 136 137 ! energy weighted density matrices in k-space 138 CALL kpoint_density_matrices(kpoints, energy_weighted=.TRUE.) 139 ! energy weighted density matrices in real space 140 CALL kpoint_density_transform(kpoints, matrix_w_kp, .TRUE., & 141 matrix_s_kp(1, 1)%matrix, sab_nl, fmwork) 142 143 DO is = 1, SIZE(fmwork) 144 CALL cp_fm_release(fmwork(is)%matrix) 145 END DO 146 DEALLOCATE (fmwork) 147 148 ELSE 149 150 NULLIFY (dft_control, rho_ao) 151 CALL get_qs_env(qs_env, & 152 matrix_w=matrix_w, & 153 matrix_ks=matrix_ks, & 154 matrix_s=matrix_s, & 155 matrix_w_mp2=matrix_w_mp2, & 156 mo_derivs=mo_derivs, & 157 scf_control=scf_control, & 158 mos=mos, & 159 rho=rho, & 160 dft_control=dft_control) 161 162 CALL qs_rho_get(rho, rho_ao=rho_ao) 163 164 nspin = SIZE(mos) 165 DO ispin = 1, nspin 166 mo_set => mos(ispin)%mo_set 167 IF (dft_control%roks) THEN 168 IF (scf_control%use_ot) THEN 169 IF (ispin > 1) THEN 170 ! not very elegant, indeed ... 171 CALL dbcsr_set(matrix_w(ispin)%matrix, 0.0_dp) 172 ELSE 173 CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, & 174 matrix_w(ispin)%matrix, matrix_s(1)%matrix) 175 END IF 176 ELSE 177 CALL calculate_w_matrix(mo_set=mo_set, & 178 matrix_ks=matrix_ks(ispin)%matrix, & 179 matrix_p=rho_ao(ispin)%matrix, & 180 matrix_w=matrix_w(ispin)%matrix) 181 END IF 182 ELSE 183 IF (scf_control%use_ot) THEN 184 CALL calculate_w_matrix_ot(mo_set, mo_derivs(ispin)%matrix, & 185 matrix_w(ispin)%matrix, matrix_s(1)%matrix) 186 ELSE 187 CALL calculate_w_matrix(mo_set, matrix_w(ispin)%matrix) 188 END IF 189 END IF 190 ! if MP2 time to update the W matrix with the MP2 contribution 191 IF (ASSOCIATED(qs_env%mp2_env)) THEN 192 CALL dbcsr_add(matrix_w(ispin)%matrix, matrix_w_mp2(ispin)%matrix, 1.0_dp, -1.0_dp) 193 END IF 194 END DO 195 196 END IF 197 198 END IF 199 200 CALL timestop(handle) 201 202 END SUBROUTINE qs_energies_compute_matrix_w 203 204! ************************************************************************************************** 205!> \brief Refactoring of qs_energies_scf. Moves computation of properties 206!> into separate subroutine 207!> \param qs_env ... 208!> \par History 209!> 05.2013 created [Florian Schiffmann] 210! ************************************************************************************************** 211 212 SUBROUTINE qs_energies_properties(qs_env) 213 TYPE(qs_environment_type), POINTER :: qs_env 214 215 CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_properties', & 216 routineP = moduleN//':'//routineN 217 218 INTEGER :: handle 219 LOGICAL :: do_et, do_et_proj 220 TYPE(atprop_type), POINTER :: atprop 221 TYPE(dft_control_type), POINTER :: dft_control 222 TYPE(pw_env_type), POINTER :: pw_env 223 TYPE(pw_p_type) :: v_hartree_rspace 224 TYPE(qs_energy_type), POINTER :: energy 225 TYPE(section_vals_type), POINTER :: input, proj_section, rest_b_section 226 227 NULLIFY (atprop, energy, pw_env) 228 CALL timeset(routineN, handle) 229 230 NULLIFY (v_hartree_rspace%pw) 231 232 ! atomic energies using Mulliken partition 233 CALL get_qs_env(qs_env, & 234 dft_control=dft_control, & 235 input=input, & 236 atprop=atprop, & 237 energy=energy, & 238 v_hartree_rspace=v_hartree_rspace%pw, & 239 pw_env=pw_env) 240 IF (atprop%energy) THEN 241 CALL qs_energies_mulliken(qs_env) 242 IF (.NOT. dft_control%qs_control%semi_empirical .AND. & 243 .NOT. dft_control%qs_control%xtb .AND. & 244 .NOT. dft_control%qs_control%dftb) THEN 245 ! Nuclear charge correction 246 CALL integrate_v_core_rspace(v_hartree_rspace, qs_env) 247 ! Kohn-Sham Functional corrections 248 END IF 249 CALL atprop_array_add(atprop%atener, atprop%ateb) 250 CALL atprop_array_add(atprop%atener, atprop%ateself) 251 CALL atprop_array_add(atprop%atener, atprop%atexc) 252 CALL atprop_array_add(atprop%atener, atprop%atecoul) 253 CALL atprop_array_add(atprop%atener, atprop%atevdw) 254 CALL atprop_array_add(atprop%atener, atprop%ategcp) 255 CALL atprop_array_add(atprop%atener, atprop%atecc) 256 CALL atprop_array_add(atprop%atener, atprop%ate1c) 257 END IF 258 259 ! ET coupling - projection-operator approach 260 NULLIFY (proj_section) 261 proj_section => & 262 section_vals_get_subs_vals(qs_env%input, "PROPERTIES%ET_COUPLING%PROJECTION") 263 CALL section_vals_get(proj_section, explicit=do_et_proj) 264 IF (do_et_proj) THEN 265 CALL calc_et_coupling_proj(qs_env) 266 END IF 267 268 ! ********** Calculate the electron transfer coupling elements******** 269 do_et = .FALSE. 270 do_et = dft_control%qs_control%et_coupling_calc 271 IF (do_et) THEN 272 qs_env%et_coupling%energy = energy%total 273 qs_env%et_coupling%keep_matrix = .TRUE. 274 qs_env%et_coupling%first_run = .TRUE. 275 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.TRUE.) 276 qs_env%et_coupling%first_run = .FALSE. 277 IF (dft_control%qs_control%ddapc_restraint) THEN 278 rest_b_section => section_vals_get_subs_vals(input, "PROPERTIES%ET_COUPLING%DDAPC_RESTRAINT_B") 279 CALL read_ddapc_section(qs_control=dft_control%qs_control, & 280 ddapc_restraint_section=rest_b_section) 281 END IF 282 CALL scf(qs_env=qs_env) 283 qs_env%et_coupling%keep_matrix = .TRUE. 284 285 CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.TRUE.) 286 CALL calc_et_coupling(qs_env) 287 END IF 288 289 !Properties 290 IF (dft_control%do_xas_calculation) THEN 291 CALL xas(qs_env, dft_control) 292 END IF 293 294 IF (dft_control%do_xas_tdp_calculation) THEN 295 CALL xas_tdp(qs_env) 296 END IF 297 298 ! Compute Linear Response properties as post-scf 299 IF (.NOT. qs_env%linres_run) THEN 300 CALL linres_calculation_low(qs_env) 301 END IF 302 303 IF (dft_control%tddfpt2_control%enabled) THEN 304 CALL tddfpt(qs_env) 305 END IF 306 307 CALL timestop(handle) 308 309 END SUBROUTINE qs_energies_properties 310 311! ************************************************************************************************** 312!> \brief Use a simple Mulliken-like energy decomposition 313!> \param qs_env ... 314!> \date 07.2011 315!> \author JHU 316!> \version 1.0 317! ************************************************************************************************** 318 SUBROUTINE qs_energies_mulliken(qs_env) 319 320 TYPE(qs_environment_type), POINTER :: qs_env 321 322 CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_mulliken', & 323 routineP = moduleN//':'//routineN 324 325 INTEGER :: ispin 326 TYPE(atprop_type), POINTER :: atprop 327 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_ks, rho_ao 328 TYPE(qs_rho_type), POINTER :: rho 329 330 NULLIFY (atprop, matrix_h, matrix_ks, rho, rho_ao) 331 CALL get_qs_env(qs_env=qs_env, matrix_ks=matrix_ks, matrix_h=matrix_h, & 332 rho=rho, atprop=atprop) 333 CALL qs_rho_get(rho, rho_ao=rho_ao) 334 335 IF (atprop%energy) THEN 336 ! E = 0.5*Tr(H*P+F*P) 337 atprop%atener = 0._dp 338 DO ispin = 1, SIZE(rho_ao) 339 CALL atom_trace(matrix_h(1)%matrix, rho_ao(ispin)%matrix, & 340 0.5_dp, atprop%atener) 341 CALL atom_trace(matrix_ks(ispin)%matrix, rho_ao(ispin)%matrix, & 342 0.5_dp, atprop%atener) 343 END DO 344 END IF 345 346 END SUBROUTINE qs_energies_mulliken 347 348! ************************************************************************************************** 349!> \brief Compute partial trace of product of two matrices 350!> \param amat ... 351!> \param bmat ... 352!> \param factor ... 353!> \param atrace ... 354!> \par History 355!> 06.2004 created [Joost VandeVondele] 356!> \note 357!> charges are computed per spin in the LSD case 358! ************************************************************************************************** 359 SUBROUTINE atom_trace(amat, bmat, factor, atrace) 360 TYPE(dbcsr_type), POINTER :: amat, bmat 361 REAL(kind=dp), INTENT(IN) :: factor 362 REAL(KIND=dp), DIMENSION(:), POINTER :: atrace 363 364 CHARACTER(len=*), PARAMETER :: routineN = 'atom_trace', routineP = moduleN//':'//routineN 365 366 INTEGER :: blk, iblock_col, iblock_row, nblock 367 LOGICAL :: found 368 REAL(kind=dp) :: btr, mult 369 REAL(KIND=dp), DIMENSION(:, :), POINTER :: a_block, b_block 370 TYPE(dbcsr_iterator_type) :: iter 371 372 CALL dbcsr_get_info(bmat, nblkrows_total=nblock) 373 CPASSERT(nblock == SIZE(atrace)) 374 375 CALL dbcsr_iterator_start(iter, bmat) 376 DO WHILE (dbcsr_iterator_blocks_left(iter)) 377 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, b_block, blk) 378 CALL dbcsr_get_block_p(matrix=amat, & 379 row=iblock_row, col=iblock_col, BLOCK=a_block, found=found) 380 381 ! we can cycle if a block is not present 382 IF (.NOT. (ASSOCIATED(b_block) .AND. ASSOCIATED(a_block))) CYCLE 383 384 IF (iblock_row .EQ. iblock_col) THEN 385 mult = 0.5_dp ! avoid double counting of diagonal blocks 386 ELSE 387 mult = 1.0_dp 388 ENDIF 389 btr = factor*mult*SUM(a_block*b_block) 390 atrace(iblock_row) = atrace(iblock_row) + btr 391 atrace(iblock_col) = atrace(iblock_col) + btr 392 393 ENDDO 394 CALL dbcsr_iterator_stop(iter) 395 396 END SUBROUTINE atom_trace 397 398! ************************************************************************************************** 399!> \brief Enters the mp2 part of cp2k 400!> \param qs_env ... 401!> \param calc_forces ... 402! ************************************************************************************************** 403 404 SUBROUTINE qs_energies_mp2(qs_env, calc_forces) 405 TYPE(qs_environment_type), POINTER :: qs_env 406 LOGICAL, INTENT(IN) :: calc_forces 407 408 CHARACTER(len=*), PARAMETER :: routineN = 'qs_energies_mp2', & 409 routineP = moduleN//':'//routineN 410 411 LOGICAL :: should_stop 412 413! Compute MP2 energy 414 415 IF (ASSOCIATED(qs_env%mp2_env)) THEN 416 417 CALL external_control(should_stop, "MP2", target_time=qs_env%target_time, & 418 start_time=qs_env%start_time) 419 420 CALL mp2_main(qs_env=qs_env, calc_forces=calc_forces) 421 ENDIF 422 423 END SUBROUTINE qs_energies_mp2 424 425END MODULE qs_energy_utils 426