1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines for a Kim-Gordon-like partitioning into molecular subunits 8!> \par History 9!> 2012.06 created [Martin Haeufel] 10!> \author Martin Haeufel and Florian Schiffmann 11! ************************************************************************************************** 12MODULE kg_correction 13 USE atomic_kind_types, ONLY: atomic_kind_type 14 USE cp_control_types, ONLY: dft_control_type 15 USE cp_para_types, ONLY: cp_para_env_type 16 USE dbcsr_api, ONLY: dbcsr_add,& 17 dbcsr_dot,& 18 dbcsr_p_type 19 USE input_constants, ONLY: kg_tnadd_atomic,& 20 kg_tnadd_embed,& 21 kg_tnadd_embed_ri,& 22 kg_tnadd_none 23 USE kg_environment_types, ONLY: kg_environment_type 24 USE kinds, ONLY: dp 25 USE lri_environment_methods, ONLY: calculate_lri_densities,& 26 lri_kg_rho_update 27 USE lri_environment_types, ONLY: lri_density_type,& 28 lri_environment_type,& 29 lri_kind_type 30 USE lri_forces, ONLY: calculate_lri_forces 31 USE lri_ks_methods, ONLY: calculate_lri_ks_matrix 32 USE message_passing, ONLY: mp_sum 33 USE pw_env_types, ONLY: pw_env_get,& 34 pw_env_type 35 USE pw_pool_types, ONLY: pw_pool_give_back_pw,& 36 pw_pool_type 37 USE pw_types, ONLY: pw_p_type 38 USE qs_environment_types, ONLY: get_qs_env,& 39 qs_environment_type 40 USE qs_integrate_potential, ONLY: integrate_v_rspace,& 41 integrate_v_rspace_one_center 42 USE qs_ks_types, ONLY: qs_ks_env_type 43 USE qs_rho_methods, ONLY: qs_rho_rebuild,& 44 qs_rho_update_rho 45 USE qs_rho_types, ONLY: qs_rho_create,& 46 qs_rho_get,& 47 qs_rho_release,& 48 qs_rho_set,& 49 qs_rho_type 50 USE qs_vxc, ONLY: qs_vxc_create 51 USE virial_types, ONLY: virial_type 52#include "./base/base_uses.f90" 53 54 IMPLICIT NONE 55 56 PRIVATE 57 58 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kg_correction' 59 60 PUBLIC :: kg_ekin_subset 61 62CONTAINS 63 64! ************************************************************************************************** 65!> \brief Calculates the subsystem Hohenberg-Kohn kinetic energy and the forces 66!> \param qs_env ... 67!> \param ks_matrix ... 68!> \param ekin_mol ... 69!> \param calc_force ... 70!> \par History 71!> 2012.06 created [Martin Haeufel] 72!> 2014.01 added atomic potential option [JGH] 73!> \author Martin Haeufel and Florian Schiffmann 74! ************************************************************************************************** 75 SUBROUTINE kg_ekin_subset(qs_env, ks_matrix, ekin_mol, calc_force) 76 TYPE(qs_environment_type), POINTER :: qs_env 77 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix 78 REAL(KIND=dp), INTENT(out) :: ekin_mol 79 LOGICAL :: calc_force 80 81 CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_ekin_subset', routineP = moduleN//':'//routineN 82 83 LOGICAL :: lrigpw 84 TYPE(dft_control_type), POINTER :: dft_control 85 TYPE(kg_environment_type), POINTER :: kg_env 86 87 CALL get_qs_env(qs_env, kg_env=kg_env, dft_control=dft_control) 88 lrigpw = dft_control%qs_control%lrigpw 89 90 IF (kg_env%tnadd_method == kg_tnadd_embed) THEN 91 IF (lrigpw) THEN 92 CALL kg_ekin_embed_lri(qs_env, kg_env, ks_matrix, ekin_mol, calc_force) 93 ELSE 94 CALL kg_ekin_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force) 95 END IF 96 ELSE IF (kg_env%tnadd_method == kg_tnadd_embed_ri) THEN 97 CALL kg_ekin_ri_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force) 98 ELSE IF (kg_env%tnadd_method == kg_tnadd_atomic) THEN 99 CALL kg_ekin_atomic(qs_env, ks_matrix, ekin_mol) 100 ELSE IF (kg_env%tnadd_method == kg_tnadd_none) THEN 101 ekin_mol = 0.0_dp 102 ELSE 103 CPABORT("Unknown KG embedding method") 104 END IF 105 106 END SUBROUTINE kg_ekin_subset 107 108! ************************************************************************************************** 109!> \brief ... 110!> \param qs_env ... 111!> \param kg_env ... 112!> \param ks_matrix ... 113!> \param ekin_mol ... 114!> \param calc_force ... 115! ************************************************************************************************** 116 SUBROUTINE kg_ekin_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force) 117 TYPE(qs_environment_type), POINTER :: qs_env 118 TYPE(kg_environment_type), POINTER :: kg_env 119 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix 120 REAL(KIND=dp), INTENT(out) :: ekin_mol 121 LOGICAL :: calc_force 122 123 CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_ekin_embed', routineP = moduleN//':'//routineN 124 125 INTEGER :: handle, ispin, isub, natom, nspins 126 LOGICAL :: use_virial 127 REAL(KIND=dp) :: ekin_imol 128 REAL(KIND=dp), DIMENSION(3, 3) :: xcvirial 129 TYPE(cp_para_env_type), POINTER :: para_env 130 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: density_matrix 131 TYPE(dft_control_type), POINTER :: dft_control 132 TYPE(pw_env_type), POINTER :: pw_env 133 TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau 134 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 135 TYPE(qs_ks_env_type), POINTER :: ks_env 136 TYPE(qs_rho_type), POINTER :: old_rho, rho_struct 137 TYPE(virial_type), POINTER :: virial 138 139 CALL timeset(routineN, handle) 140 141 NULLIFY (vxc_rho, vxc_tau, old_rho, rho_struct, ks_env) 142 143 ! get set of molecules, natom, dft_control, pw_env 144 CALL get_qs_env(qs_env, & 145 ks_env=ks_env, & 146 rho=old_rho, & 147 natom=natom, & 148 dft_control=dft_control, & 149 virial=virial, & 150 para_env=para_env, & 151 pw_env=pw_env) 152 153 nspins = dft_control%nspins 154 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 155 use_virial = use_virial .AND. calc_force 156 157 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 158 159 ! get the density matrix 160 CALL qs_rho_get(old_rho, rho_ao=density_matrix) 161 ! allocate and initialize the density 162 CALL qs_rho_create(rho_struct) 163 ! set the density matrix to the blocked matrix 164 CALL qs_rho_set(rho_struct, rho_ao=density_matrix) ! blocked_matrix 165 CALL qs_rho_rebuild(rho_struct, qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.) 166 167 ! full density kinetic energy term 168 CALL qs_rho_update_rho(rho_struct, qs_env) 169 ekin_imol = 0.0_dp 170 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, & 171 vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol) 172 IF (ASSOCIATED(vxc_tau)) THEN 173 CPABORT(" KG with meta-kinetic energy functionals not implemented") 174 END IF 175 DO ispin = 1, nspins 176 vxc_rho(ispin)%pw%cr3d = vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol 177 CALL integrate_v_rspace(v_rspace=vxc_rho(ispin), & 178 pmat=density_matrix(ispin), hmat=ks_matrix(ispin), & 179 qs_env=qs_env, calculate_forces=calc_force) 180 CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw) 181 END DO 182 DEALLOCATE (vxc_rho) 183 ekin_mol = -ekin_imol 184 xcvirial(1:3, 1:3) = 0.0_dp 185 IF (use_virial) xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) - virial%pv_xc(1:3, 1:3) 186 187 ! loop over all subsets 188 DO isub = 1, kg_env%nsubsets 189 ! calculate the densities for the given blocked density matrix - pass the subset task_list 190 CALL qs_rho_update_rho(rho_struct, qs_env, & 191 task_list_external=kg_env%subset(isub)%task_list) 192 193 ekin_imol = 0.0_dp 194 ! calc Hohenberg-Kohn kin. energy of the density corresp. to the remaining molecular block(s) 195 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, & 196 vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol) 197 ekin_mol = ekin_mol + ekin_imol 198 199 DO ispin = 1, nspins 200 vxc_rho(ispin)%pw%cr3d = -vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol 201 CALL integrate_v_rspace(v_rspace=vxc_rho(ispin), & 202 pmat=density_matrix(ispin), & 203 hmat=ks_matrix(ispin), & 204 qs_env=qs_env, & 205 calculate_forces=calc_force, & 206 task_list_external=kg_env%subset(isub)%task_list) 207 ! clean up vxc_rho 208 CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw) 209 END DO 210 DEALLOCATE (vxc_rho) 211 212 IF (use_virial) THEN 213 xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) + virial%pv_xc(1:3, 1:3) 214 END IF 215 216 END DO 217 218 IF (use_virial) THEN 219 virial%pv_xc(1:3, 1:3) = xcvirial(1:3, 1:3) 220 END IF 221 222 ! clean up rho_struct 223 CALL qs_rho_set(rho_struct, rho_ao=Null()) 224 CALL qs_rho_release(rho_struct) 225 226 CALL timestop(handle) 227 228 END SUBROUTINE kg_ekin_embed 229 230! ************************************************************************************************** 231!> \brief ... 232!> \param qs_env ... 233!> \param kg_env ... 234!> \param ks_matrix ... 235!> \param ekin_mol ... 236!> \param calc_force ... 237! ************************************************************************************************** 238 SUBROUTINE kg_ekin_embed_lri(qs_env, kg_env, ks_matrix, ekin_mol, calc_force) 239 TYPE(qs_environment_type), POINTER :: qs_env 240 TYPE(kg_environment_type), POINTER :: kg_env 241 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix 242 REAL(KIND=dp), INTENT(out) :: ekin_mol 243 LOGICAL :: calc_force 244 245 CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_ekin_embed_lri', & 246 routineP = moduleN//':'//routineN 247 248 INTEGER :: color, handle, iatom, ikind, imol, & 249 ispin, isub, natom, nkind, nspins 250 INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist 251 LOGICAL :: use_virial 252 REAL(KIND=dp) :: ekin_imol 253 REAL(KIND=dp), DIMENSION(3, 3) :: xcvirial 254 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 255 TYPE(cp_para_env_type), POINTER :: para_env 256 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: density_matrix, ksmat 257 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: pmat 258 TYPE(dft_control_type), POINTER :: dft_control 259 TYPE(lri_density_type), POINTER :: lri_density 260 TYPE(lri_environment_type), POINTER :: lri_env 261 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int 262 TYPE(pw_env_type), POINTER :: pw_env 263 TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau 264 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 265 TYPE(qs_ks_env_type), POINTER :: ks_env 266 TYPE(qs_rho_type), POINTER :: old_rho, rho_struct 267 TYPE(virial_type), POINTER :: virial 268 269 CALL timeset(routineN, handle) 270 271 NULLIFY (vxc_rho, vxc_tau, old_rho, rho_struct, ks_env) 272 273 CALL get_qs_env(qs_env, dft_control=dft_control) 274 275 ! get set of molecules, natom, dft_control, pw_env 276 CALL get_qs_env(qs_env, & 277 ks_env=ks_env, & 278 rho=old_rho, & 279 natom=natom, & 280 dft_control=dft_control, & 281 virial=virial, & 282 para_env=para_env, & 283 pw_env=pw_env) 284 285 nspins = dft_control%nspins 286 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 287 use_virial = use_virial .AND. calc_force 288 289 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 290 291 ! get the density matrix 292 CALL qs_rho_get(old_rho, rho_ao=density_matrix) 293 ! allocate and initialize the density 294 CALL qs_rho_create(rho_struct) 295 ! set the density matrix to the blocked matrix 296 CALL qs_rho_set(rho_struct, rho_ao=density_matrix) ! blocked_matrix 297 CALL qs_rho_rebuild(rho_struct, qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.) 298 299 CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density, nkind=nkind) 300 IF (lri_env%exact_1c_terms) THEN 301 CPABORT(" KG with LRI and exact one-center terms not implemented") 302 END IF 303 ALLOCATE (atomlist(natom)) 304 DO ispin = 1, nspins 305 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 306 DO ikind = 1, nkind 307 lri_v_int(ikind)%v_int = 0.0_dp 308 IF (calc_force) THEN 309 lri_v_int(ikind)%v_dadr = 0.0_dp 310 lri_v_int(ikind)%v_dfdr = 0.0_dp 311 END IF 312 END DO 313 END DO 314 315 ! full density kinetic energy term 316 atomlist = 1 317 CALL lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist) 318 ekin_imol = 0.0_dp 319 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, & 320 vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol) 321 IF (ASSOCIATED(vxc_tau)) THEN 322 CPABORT(" KG with meta-kinetic energy functionals not implemented") 323 END IF 324 DO ispin = 1, nspins 325 vxc_rho(ispin)%pw%cr3d = vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol 326 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 327 CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, lri_v_int, calc_force, "LRI_AUX") 328 CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw) 329 END DO 330 DEALLOCATE (vxc_rho) 331 ekin_mol = -ekin_imol 332 xcvirial(1:3, 1:3) = 0.0_dp 333 IF (use_virial) xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) - virial%pv_xc(1:3, 1:3) 334 335 ! loop over all subsets 336 DO isub = 1, kg_env%nsubsets 337 atomlist = 0 338 DO iatom = 1, natom 339 imol = kg_env%atom_to_molecule(iatom) 340 color = kg_env%subset_of_mol(imol) 341 IF (color == isub) atomlist(iatom) = 1 342 END DO 343 CALL lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist) 344 345 ekin_imol = 0.0_dp 346 ! calc Hohenberg-Kohn kin. energy of the density corresp. to the remaining molecular block(s) 347 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, & 348 vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol) 349 ekin_mol = ekin_mol + ekin_imol 350 351 DO ispin = 1, nspins 352 vxc_rho(ispin)%pw%cr3d = -vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol 353 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 354 CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, & 355 lri_v_int, calc_force, & 356 "LRI_AUX", atomlist=atomlist) 357 ! clean up vxc_rho 358 CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw) 359 END DO 360 DEALLOCATE (vxc_rho) 361 362 IF (use_virial) THEN 363 xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) + virial%pv_xc(1:3, 1:3) 364 END IF 365 366 END DO 367 368 IF (use_virial) THEN 369 virial%pv_xc(1:3, 1:3) = xcvirial(1:3, 1:3) 370 END IF 371 372 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set) 373 ALLOCATE (ksmat(1)) 374 DO ispin = 1, nspins 375 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 376 DO ikind = 1, nkind 377 CALL mp_sum(lri_v_int(ikind)%v_int, para_env%group) 378 END DO 379 ksmat(1)%matrix => ks_matrix(ispin)%matrix 380 CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set) 381 END DO 382 IF (calc_force) THEN 383 pmat(1:nspins, 1:1) => density_matrix(1:nspins) 384 CALL calculate_lri_forces(lri_env, lri_density, qs_env, pmat, atomic_kind_set) 385 ENDIF 386 DEALLOCATE (atomlist, ksmat) 387 388 ! clean up rho_struct 389 CALL qs_rho_set(rho_struct, rho_ao=Null()) 390 CALL qs_rho_release(rho_struct) 391 392 CALL timestop(handle) 393 394 END SUBROUTINE kg_ekin_embed_lri 395 396! ************************************************************************************************** 397!> \brief ... 398!> \param qs_env ... 399!> \param kg_env ... 400!> \param ks_matrix ... 401!> \param ekin_mol ... 402!> \param calc_force ... 403! ************************************************************************************************** 404 SUBROUTINE kg_ekin_ri_embed(qs_env, kg_env, ks_matrix, ekin_mol, calc_force) 405 TYPE(qs_environment_type), POINTER :: qs_env 406 TYPE(kg_environment_type), POINTER :: kg_env 407 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix 408 REAL(KIND=dp), INTENT(out) :: ekin_mol 409 LOGICAL :: calc_force 410 411 CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_ekin_ri_embed', & 412 routineP = moduleN//':'//routineN 413 414 INTEGER :: color, handle, iatom, ikind, imol, & 415 ispin, isub, natom, nkind, nspins 416 INTEGER, ALLOCATABLE, DIMENSION(:) :: atomlist 417 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 418 LOGICAL :: use_virial 419 REAL(KIND=dp) :: ekin_imol 420 REAL(KIND=dp), DIMENSION(3, 3) :: xcvirial 421 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 422 TYPE(cp_para_env_type), POINTER :: para_env 423 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat 424 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: density_matrix 425 TYPE(dft_control_type), POINTER :: dft_control 426 TYPE(lri_density_type), POINTER :: lri_density 427 TYPE(lri_environment_type), POINTER :: lri_env 428 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int 429 TYPE(pw_env_type), POINTER :: pw_env 430 TYPE(pw_p_type), DIMENSION(:), POINTER :: vxc_rho, vxc_tau 431 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 432 TYPE(qs_ks_env_type), POINTER :: ks_env 433 TYPE(qs_rho_type), POINTER :: rho, rho_struct 434 TYPE(virial_type), POINTER :: virial 435 436 CALL timeset(routineN, handle) 437 438 CALL get_qs_env(qs_env, & 439 ks_env=ks_env, & 440 rho=rho, & 441 natom=natom, & 442 nkind=nkind, & 443 dft_control=dft_control, & 444 virial=virial, & 445 para_env=para_env, & 446 pw_env=pw_env) 447 448 nspins = dft_control%nspins 449 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) 450 use_virial = use_virial .AND. calc_force 451 452 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 453 454 ! get the density matrix 455 CALL qs_rho_get(rho, rho_ao_kp=density_matrix) 456 ! allocate and initialize the density 457 NULLIFY (rho_struct) 458 CALL qs_rho_create(rho_struct) 459 ! set the density matrix to the blocked matrix 460 CALL qs_rho_set(rho_struct, rho_ao_kp=density_matrix) 461 CALL qs_rho_rebuild(rho_struct, qs_env, rebuild_ao=.FALSE., rebuild_grids=.TRUE.) 462 463 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set) 464 ALLOCATE (cell_to_index(1, 1, 1)) 465 cell_to_index(1, 1, 1) = 1 466 lri_env => kg_env%lri_env 467 lri_density => kg_env%lri_density 468 CALL calculate_lri_densities(lri_env, lri_density, qs_env, density_matrix, cell_to_index, & 469 rho_struct, atomic_kind_set, para_env) 470 kg_env%lri_density => lri_density 471 ! full density kinetic energy term 472 ekin_imol = 0.0_dp 473 NULLIFY (vxc_rho, vxc_tau) 474 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, & 475 vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol) 476 IF (ASSOCIATED(vxc_tau)) THEN 477 CPABORT(" KG with meta-kinetic energy functionals not implemented") 478 END IF 479 DO ispin = 1, nspins 480 vxc_rho(ispin)%pw%cr3d = vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol 481 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 482 CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, lri_v_int, calc_force, "LRI_AUX") 483 CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw) 484 END DO 485 DEALLOCATE (vxc_rho) 486 ekin_mol = -ekin_imol 487 xcvirial(1:3, 1:3) = 0.0_dp 488 IF (use_virial) xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) - virial%pv_xc(1:3, 1:3) 489!deb 490! WRITE(6,*) " E KIN (full) ",-ekin_mol 491!deb 492 493 ! loop over all subsets 494 ALLOCATE (atomlist(natom)) 495 DO isub = 1, kg_env%nsubsets 496 atomlist = 0 497 DO iatom = 1, natom 498 imol = kg_env%atom_to_molecule(iatom) 499 color = kg_env%subset_of_mol(imol) 500 IF (color == isub) atomlist(iatom) = 1 501 END DO 502 CALL lri_kg_rho_update(rho_struct, qs_env, lri_env, lri_density, atomlist) 503 504 ekin_imol = 0.0_dp 505 ! calc Hohenberg-Kohn kin. energy of the density corresp. to the remaining molecular block(s) 506 NULLIFY (vxc_rho, vxc_tau) 507 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=kg_env%xc_section_kg, & 508 vxc_rho=vxc_rho, vxc_tau=vxc_tau, exc=ekin_imol) 509 ekin_mol = ekin_mol + ekin_imol 510!deb 511! WRITE(6,*) " E KIN (molecule) ",isub,ekin_imol 512!deb 513 514 DO ispin = 1, nspins 515 vxc_rho(ispin)%pw%cr3d = -vxc_rho(ispin)%pw%cr3d*vxc_rho(ispin)%pw%pw_grid%dvol 516 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 517 CALL integrate_v_rspace_one_center(vxc_rho(ispin), qs_env, & 518 lri_v_int, calc_force, & 519 "LRI_AUX", atomlist=atomlist) 520 ! clean up vxc_rho 521 CALL pw_pool_give_back_pw(auxbas_pw_pool, vxc_rho(ispin)%pw) 522 END DO 523 DEALLOCATE (vxc_rho) 524 525 IF (use_virial) THEN 526 xcvirial(1:3, 1:3) = xcvirial(1:3, 1:3) + virial%pv_xc(1:3, 1:3) 527 END IF 528 529 END DO 530 531 IF (use_virial) THEN 532 virial%pv_xc(1:3, 1:3) = xcvirial(1:3, 1:3) 533 END IF 534 535 ALLOCATE (ksmat(1)) 536 DO ispin = 1, nspins 537 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds 538 DO ikind = 1, nkind 539 CALL mp_sum(lri_v_int(ikind)%v_int, para_env%group) 540 END DO 541 ksmat(1)%matrix => ks_matrix(ispin)%matrix 542 CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set) 543 END DO 544 IF (calc_force) THEN 545 CALL calculate_lri_forces(lri_env, lri_density, qs_env, density_matrix, atomic_kind_set) 546 ENDIF 547 DEALLOCATE (atomlist, ksmat) 548 549 ! clean up rho_struct 550 CALL qs_rho_set(rho_struct, rho_ao=Null()) 551 CALL qs_rho_release(rho_struct) 552 DEALLOCATE (cell_to_index) 553 554 CALL timestop(handle) 555 556 END SUBROUTINE kg_ekin_ri_embed 557 558! ************************************************************************************************** 559!> \brief ... 560!> \param qs_env ... 561!> \param ks_matrix ... 562!> \param ekin_mol ... 563! ************************************************************************************************** 564 SUBROUTINE kg_ekin_atomic(qs_env, ks_matrix, ekin_mol) 565 TYPE(qs_environment_type), POINTER :: qs_env 566 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix 567 REAL(KIND=dp), INTENT(out) :: ekin_mol 568 569 CHARACTER(LEN=*), PARAMETER :: routineN = 'kg_ekin_atomic', routineP = moduleN//':'//routineN 570 571 INTEGER :: handle, ispin, nspins 572 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: density_matrix, tnadd_matrix 573 TYPE(kg_environment_type), POINTER :: kg_env 574 TYPE(qs_rho_type), POINTER :: rho 575 576 NULLIFY (rho, kg_env, density_matrix, tnadd_matrix) 577 578 CALL timeset(routineN, handle) 579 CALL get_qs_env(qs_env, kg_env=kg_env, rho=rho) 580 581 nspins = SIZE(ks_matrix) 582 ! get the density matrix 583 CALL qs_rho_get(rho, rho_ao=density_matrix) 584 ! get the tnadd matrix 585 tnadd_matrix => kg_env%tnadd_mat 586 587 ekin_mol = 0.0_dp 588 DO ispin = 1, nspins 589 CALL dbcsr_dot(tnadd_matrix(1)%matrix, density_matrix(ispin)%matrix, ekin_mol) 590 CALL dbcsr_add(ks_matrix(ispin)%matrix, tnadd_matrix(1)%matrix, & 591 alpha_scalar=1.0_dp, beta_scalar=1.0_dp) 592 END DO 593 ! definition is inverted (see qs_ks_methods) 594 ekin_mol = -ekin_mol 595 596 CALL timestop(handle) 597 598 END SUBROUTINE kg_ekin_atomic 599 600END MODULE kg_correction 601