1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief methods of the rho structure (defined in qs_rho_types) 8!> \par History 9!> 08.2002 created [fawzi] 10!> 08.2014 kpoints [JGH] 11!> \author Fawzi Mohamed 12! ************************************************************************************************** 13MODULE qs_rho_methods 14 USE atomic_kind_types, ONLY: atomic_kind_type 15 USE cp_control_types, ONLY: dft_control_type 16 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 17 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,& 18 dbcsr_deallocate_matrix_set 19 USE cp_log_handling, ONLY: cp_to_string 20 USE cp_para_types, ONLY: cp_para_env_type 21 USE dbcsr_api, ONLY: dbcsr_copy,& 22 dbcsr_create,& 23 dbcsr_p_type,& 24 dbcsr_set,& 25 dbcsr_type,& 26 dbcsr_type_symmetric 27 USE kinds, ONLY: default_string_length,& 28 dp 29 USE kpoint_types, ONLY: get_kpoint_info,& 30 kpoint_type 31 USE lri_environment_methods, ONLY: calculate_lri_densities 32 USE lri_environment_types, ONLY: lri_density_type,& 33 lri_environment_type 34 USE pw_env_types, ONLY: pw_env_get,& 35 pw_env_type 36 USE pw_methods, ONLY: pw_zero 37 USE pw_pool_types, ONLY: pw_pool_create_pw,& 38 pw_pool_type 39 USE pw_types, ONLY: COMPLEXDATA1D,& 40 REALDATA3D,& 41 REALSPACE,& 42 RECIPROCALSPACE,& 43 pw_p_type,& 44 pw_release 45 USE qs_collocate_density, ONLY: calculate_drho_elec,& 46 calculate_rho_elec 47 USE qs_environment_types, ONLY: get_qs_env,& 48 qs_environment_type,& 49 set_qs_env 50 USE qs_ks_types, ONLY: get_ks_env,& 51 qs_ks_env_type 52 USE qs_local_rho_types, ONLY: local_rho_type 53 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type 54 USE qs_rho_atom_methods, ONLY: calculate_rho_atom_coeff 55 USE qs_rho_types, ONLY: qs_rho_clear,& 56 qs_rho_get,& 57 qs_rho_set,& 58 qs_rho_type 59 USE ri_environment_methods, ONLY: calculate_ri_densities 60 USE task_list_types, ONLY: task_list_type 61#include "./base/base_uses.f90" 62 63 IMPLICIT NONE 64 PRIVATE 65 66 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE. 67 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho_methods' 68 69 PUBLIC :: qs_rho_update_rho, qs_rho_rebuild, duplicate_rho_type 70 71CONTAINS 72 73! ************************************************************************************************** 74!> \brief rebuilds rho (if necessary allocating and initializing it) 75!> \param rho the rho type to rebuild (defaults to qs_env%rho) 76!> \param qs_env the environment to which rho belongs 77!> \param rebuild_ao if it is necessary to rebuild rho_ao. Defaults to true. 78!> \param rebuild_grids if it in necessary to rebuild rho_r and rho_g. 79!> Defaults to false. 80!> \param admm (use aux_fit basis) 81!> \param pw_env_external external plane wave environment 82!> \par History 83!> 11.2002 created replacing qs_rho_create and qs_env_rebuild_rho[fawzi] 84!> \author Fawzi Mohamed 85!> \note 86!> needs updated pw pools, s, s_mstruct and h in qs_env. 87!> The use of p to keep the structure of h (needed for the forces) 88!> is ugly and should be removed. 89!> Change so that it does not allocate a subcomponent if it is not 90!> associated and not requested? 91! ************************************************************************************************** 92 SUBROUTINE qs_rho_rebuild(rho, qs_env, rebuild_ao, rebuild_grids, admm, pw_env_external) 93 TYPE(qs_rho_type), POINTER :: rho 94 TYPE(qs_environment_type), POINTER :: qs_env 95 LOGICAL, INTENT(in), OPTIONAL :: rebuild_ao, rebuild_grids, admm 96 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external 97 98 CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_rho_rebuild', routineP = moduleN//':'//routineN 99 100 CHARACTER(LEN=default_string_length) :: headline 101 INTEGER :: handle, i, ic, nimg, nspins 102 LOGICAL :: do_kpoints, my_admm, my_rebuild_ao, & 103 my_rebuild_grids 104 REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r 105 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 106 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, rho_ao_kp 107 TYPE(dbcsr_type), POINTER :: refmatrix, tmatrix 108 TYPE(dft_control_type), POINTER :: dft_control 109 TYPE(kpoint_type), POINTER :: kpoints 110 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 111 POINTER :: sab_orb 112 TYPE(pw_env_type), POINTER :: pw_env 113 TYPE(pw_p_type), DIMENSION(:), POINTER :: drho_g, drho_r, rho_g, rho_r, tau_g, & 114 tau_r 115 TYPE(pw_p_type), POINTER :: rho_r_sccs 116 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 117 118 CALL timeset(routineN, handle) 119 120 NULLIFY (pw_env, auxbas_pw_pool, matrix_s, matrix_s_kp, dft_control) 121 NULLIFY (tot_rho_r, rho_ao_kp, rho_r, rho_g, drho_r, drho_g, tau_r, tau_g) 122 NULLIFY (rho_r_sccs) 123 NULLIFY (sab_orb) 124 my_rebuild_ao = .TRUE. 125 my_rebuild_grids = .TRUE. 126 my_admm = .FALSE. 127 IF (PRESENT(rebuild_ao)) my_rebuild_ao = rebuild_ao 128 IF (PRESENT(rebuild_grids)) my_rebuild_grids = rebuild_grids 129 IF (PRESENT(admm)) my_admm = admm 130 131 CALL get_qs_env(qs_env, & 132 kpoints=kpoints, & 133 do_kpoints=do_kpoints, & 134 pw_env=pw_env, & 135 dft_control=dft_control) 136 IF (PRESENT(pw_env_external)) & 137 pw_env => pw_env_external 138 139 nimg = dft_control%nimages 140 141 IF (my_admm) THEN 142 CPASSERT(.NOT. do_kpoints) 143 CALL get_qs_env(qs_env, matrix_s_aux_fit=matrix_s, sab_aux_fit=sab_orb) 144 refmatrix => matrix_s(1)%matrix 145 ELSE 146 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s_kp, sab_orb=sab_orb) 147 refmatrix => matrix_s_kp(1, 1)%matrix 148 END IF 149 150 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 151 nspins = dft_control%nspins 152 153 IF (.NOT. ASSOCIATED(rho)) CPABORT("rho not associated") 154 CALL qs_rho_get(rho, & 155 tot_rho_r=tot_rho_r, & 156 rho_ao_kp=rho_ao_kp, & 157 rho_r=rho_r, & 158 rho_g=rho_g, & 159 drho_r=drho_r, & 160 drho_g=drho_g, & 161 tau_r=tau_r, & 162 tau_g=tau_g, & 163 rho_r_sccs=rho_r_sccs) 164 165 IF (.NOT. ASSOCIATED(tot_rho_r)) THEN 166 ALLOCATE (tot_rho_r(nspins)) 167 tot_rho_r = 0.0_dp 168 CALL qs_rho_set(rho, tot_rho_r=tot_rho_r) 169 END IF 170 171 ! rho_ao 172 IF (my_rebuild_ao .OR. (.NOT. ASSOCIATED(rho_ao_kp))) THEN 173 IF (ASSOCIATED(rho_ao_kp)) & 174 CALL dbcsr_deallocate_matrix_set(rho_ao_kp) 175 ! Create a new density matrix set 176 CALL dbcsr_allocate_matrix_set(rho_ao_kp, nspins, nimg) 177 CALL qs_rho_set(rho, rho_ao_kp=rho_ao_kp) 178 DO i = 1, nspins 179 DO ic = 1, nimg 180 IF (nspins > 1) THEN 181 IF (i == 1) THEN 182 headline = "DENSITY MATRIX FOR ALPHA SPIN" 183 ELSE 184 headline = "DENSITY MATRIX FOR BETA SPIN" 185 END IF 186 ELSE 187 headline = "DENSITY MATRIX" 188 END IF 189 ALLOCATE (rho_ao_kp(i, ic)%matrix) 190 tmatrix => rho_ao_kp(i, ic)%matrix 191 CALL dbcsr_create(matrix=tmatrix, template=refmatrix, name=TRIM(headline), & 192 matrix_type=dbcsr_type_symmetric, nze=0) 193 CALL cp_dbcsr_alloc_block_from_nbl(tmatrix, sab_orb) 194 CALL dbcsr_set(tmatrix, 0.0_dp) 195 END DO 196 END DO 197 END IF 198 199 ! rho_r 200 IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(rho_r)) THEN 201 IF (ASSOCIATED(rho_r)) THEN 202 DO i = 1, SIZE(rho_r) 203 CALL pw_release(rho_r(i)%pw) 204 END DO 205 DEALLOCATE (rho_r) 206 END IF 207 ALLOCATE (rho_r(nspins)) 208 CALL qs_rho_set(rho, rho_r=rho_r) 209 DO i = 1, nspins 210 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r(i)%pw, & 211 use_data=REALDATA3D, in_space=REALSPACE) 212 END DO 213 END IF 214 215 ! rho_g 216 IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(rho_g)) THEN 217 IF (ASSOCIATED(rho_g)) THEN 218 DO i = 1, SIZE(rho_g) 219 CALL pw_release(rho_g(i)%pw) 220 END DO 221 DEALLOCATE (rho_g) 222 END IF 223 ALLOCATE (rho_g(nspins)) 224 CALL qs_rho_set(rho, rho_g=rho_g) 225 DO i = 1, nspins 226 CALL pw_pool_create_pw(auxbas_pw_pool, rho_g(i)%pw, & 227 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 228 END DO 229 END IF 230 231 ! SCCS 232 IF (dft_control%do_sccs) THEN 233 IF (my_rebuild_grids .OR. (.NOT. ASSOCIATED(rho_r_sccs))) THEN 234 IF (ASSOCIATED(rho_r_sccs)) THEN 235 CALL pw_release(rho_r_sccs%pw) 236 DEALLOCATE (rho_r_sccs) 237 END IF 238 ALLOCATE (rho_r_sccs) 239 CALL qs_rho_set(rho, rho_r_sccs=rho_r_sccs) 240 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r_sccs%pw, & 241 use_data=REALDATA3D, & 242 in_space=REALSPACE) 243 CALL pw_zero(rho_r_sccs%pw) 244 END IF 245 END IF 246 247 ! allocate drho_r and drho_g if xc_deriv_collocate 248 IF (dft_control%drho_by_collocation) THEN 249 ! drho_r 250 IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(drho_r)) THEN 251 IF (ASSOCIATED(drho_r)) THEN 252 DO i = 1, SIZE(drho_r) 253 CALL pw_release(drho_r(i)%pw) 254 END DO 255 DEALLOCATE (drho_r) 256 END IF 257 ALLOCATE (drho_r(3*nspins)) 258 CALL qs_rho_set(rho, drho_r=drho_r) 259 DO i = 1, 3*nspins 260 CALL pw_pool_create_pw(auxbas_pw_pool, drho_r(i)%pw, & 261 use_data=REALDATA3D, in_space=REALSPACE) 262 END DO 263 END IF 264 ! drho_g 265 IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(drho_g)) THEN 266 IF (ASSOCIATED(drho_g)) THEN 267 DO i = 1, SIZE(drho_g) 268 CALL pw_release(drho_g(i)%pw) 269 END DO 270 DEALLOCATE (drho_g) 271 END IF 272 ALLOCATE (drho_g(3*nspins)) 273 CALL qs_rho_set(rho, drho_g=drho_g) 274 DO i = 1, 3*nspins 275 CALL pw_pool_create_pw(auxbas_pw_pool, drho_g(i)%pw, & 276 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 277 END DO 278 END IF 279 END IF 280 281 ! allocate tau_r and tau_g if use_kinetic_energy_density 282 IF (dft_control%use_kinetic_energy_density) THEN 283 ! tau_r 284 IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(tau_r)) THEN 285 IF (ASSOCIATED(tau_r)) THEN 286 DO i = 1, SIZE(tau_r) 287 CALL pw_release(tau_r(i)%pw) 288 END DO 289 DEALLOCATE (tau_r) 290 END IF 291 ALLOCATE (tau_r(nspins)) 292 CALL qs_rho_set(rho, tau_r=tau_r) 293 DO i = 1, nspins 294 CALL pw_pool_create_pw(auxbas_pw_pool, tau_r(i)%pw, & 295 use_data=REALDATA3D, in_space=REALSPACE) 296 END DO 297 END IF 298 299 ! tau_g 300 IF (my_rebuild_grids .OR. .NOT. ASSOCIATED(tau_g)) THEN 301 IF (ASSOCIATED(tau_g)) THEN 302 DO i = 1, SIZE(tau_g) 303 CALL pw_release(tau_g(i)%pw) 304 END DO 305 DEALLOCATE (tau_g) 306 END IF 307 ALLOCATE (tau_g(nspins)) 308 CALL qs_rho_set(rho, tau_g=tau_g) 309 DO i = 1, nspins 310 CALL pw_pool_create_pw(auxbas_pw_pool, tau_g(i)%pw, & 311 use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 312 END DO 313 END IF 314 END IF ! use_kinetic_energy_density 315 316 CALL timestop(handle) 317 318 END SUBROUTINE qs_rho_rebuild 319 320! ************************************************************************************************** 321!> \brief updates rho_r and rho_g to the rho%rho_ao. 322!> if use_kinetic_energy_density also computes tau_r and tau_g 323!> \param rho_struct the rho structure that should be updated 324!> \param qs_env the qs_env rho_struct refers to 325!> the integrated charge in r space 326!> \param local_rho_set ... 327!> \param pw_env_external external plane wave environment 328!> \param task_list_external external task list 329!> \par History 330!> 08.2002 created [fawzi] 331!> \author Fawzi Mohamed 332! ************************************************************************************************** 333 SUBROUTINE qs_rho_update_rho(rho_struct, qs_env, local_rho_set, pw_env_external, task_list_external) 334 TYPE(qs_rho_type), POINTER :: rho_struct 335 TYPE(qs_environment_type), POINTER :: qs_env 336 TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set 337 TYPE(pw_env_type), OPTIONAL, POINTER :: pw_env_external 338 TYPE(task_list_type), OPTIONAL, POINTER :: task_list_external 339 340 CHARACTER(len=*), PARAMETER :: routineN = 'qs_rho_update_rho', & 341 routineP = moduleN//':'//routineN 342 343 INTEGER :: handle, img, ispin, nimg, nspins 344 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index 345 LOGICAL :: gapw, gapw_xc 346 REAL(KIND=dp) :: dum 347 REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r, tot_rho_r_xc 348 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 349 TYPE(cp_para_env_type), POINTER :: para_env 350 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao 351 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp, rho_xc_ao 352 TYPE(dft_control_type), POINTER :: dft_control 353 TYPE(kpoint_type), POINTER :: kpoints 354 TYPE(lri_density_type), POINTER :: lri_density 355 TYPE(lri_environment_type), POINTER :: lri_env 356 TYPE(pw_env_type), POINTER :: pw_env 357 TYPE(pw_p_type), DIMENSION(:), POINTER :: drho_g, drho_r, drho_xc_g, rho_g, rho_r, & 358 rho_xc_g, rho_xc_r, tau_g, tau_r, & 359 tau_xc_g, tau_xc_r 360 TYPE(pw_p_type), POINTER :: rho_r_sccs 361 TYPE(qs_ks_env_type), POINTER :: ks_env 362 TYPE(qs_rho_type), POINTER :: rho_xc 363 TYPE(task_list_type), POINTER :: task_list 364 365 CALL timeset(routineN, handle) 366 367 NULLIFY (dft_control, rho_xc, ks_env, rho_ao, rho_r, rho_g, drho_r, drho_g, tau_r, tau_g) 368 NULLIFY (rho_xc_ao, rho_xc_g, rho_xc_r, drho_xc_g, tau_xc_r, tau_xc_g, tot_rho_r, tot_rho_r_xc) 369 NULLIFY (lri_env, para_env, pw_env, atomic_kind_set) 370 371 CPASSERT(ASSOCIATED(rho_struct)) 372 373 CALL get_qs_env(qs_env, & 374 ks_env=ks_env, & 375 dft_control=dft_control, & 376 task_list=task_list, & 377 lri_env=lri_env, & 378 atomic_kind_set=atomic_kind_set, & 379 para_env=para_env, & 380 pw_env=pw_env) 381 382 CALL qs_rho_get(rho_struct, & 383 rho_r=rho_r, & 384 rho_g=rho_g, & 385 tot_rho_r=tot_rho_r, & 386 drho_r=drho_r, & 387 drho_g=drho_g, & 388 tau_r=tau_r, & 389 tau_g=tau_g, & 390 rho_r_sccs=rho_r_sccs) 391 392 IF (PRESENT(pw_env_external)) pw_env => pw_env_external 393 IF (PRESENT(task_list_external)) task_list => task_list_external 394 395 nspins = dft_control%nspins 396 nimg = dft_control%nimages 397 gapw = dft_control%qs_control%gapw 398 gapw_xc = dft_control%qs_control%gapw_xc 399 400 IF (dft_control%qs_control%semi_empirical) THEN 401 ! 402 CALL qs_rho_set(rho_struct, rho_r_valid=.FALSE., rho_g_valid=.FALSE.) 403 ELSEIF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN 404 ! 405 CALL qs_rho_set(rho_struct, rho_r_valid=.FALSE., rho_g_valid=.FALSE.) 406 ELSEIF (dft_control%qs_control%lrigpw) THEN 407 CPASSERT(.NOT. dft_control%use_kinetic_energy_density) 408 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints) 409 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index) 410 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp) 411 CALL get_qs_env(qs_env, lri_density=lri_density) 412 CALL calculate_lri_densities(lri_env, lri_density, qs_env, rho_ao_kp, cell_to_index, & 413 lri_rho_struct=rho_struct, & 414 atomic_kind_set=atomic_kind_set, & 415 para_env=para_env) 416 CALL set_qs_env(qs_env, lri_density=lri_density) 417 CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.) 418 ELSEIF (dft_control%qs_control%rigpw) THEN 419 CPASSERT(.NOT. dft_control%use_kinetic_energy_density) 420 CALL qs_rho_get(rho_struct, rho_ao=rho_ao) 421 CALL calculate_ri_densities(lri_env, qs_env, rho_ao, & 422 lri_rho_struct=rho_struct, & 423 atomic_kind_set=atomic_kind_set, & 424 para_env=para_env) 425 CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.) 426 ELSE 427 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp) 428 DO ispin = 1, nspins 429 rho_ao => rho_ao_kp(ispin, :) 430 CALL calculate_rho_elec(matrix_p_kp=rho_ao, & 431 rho=rho_r(ispin), & 432 rho_gspace=rho_g(ispin), & 433 total_rho=tot_rho_r(ispin), & 434 ks_env=ks_env, soft_valid=gapw, & 435 task_list_external=task_list, & 436 pw_env_external=pw_env) 437 END DO 438 CALL qs_rho_set(rho_struct, rho_r_valid=.TRUE., rho_g_valid=.TRUE.) 439 440 ! if needed compute also the gradient of the density 441 IF (dft_control%drho_by_collocation) THEN 442 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp) 443 CPASSERT(.NOT. PRESENT(task_list_external)) 444 CPASSERT(.NOT. PRESENT(pw_env_external)) 445 DO ispin = 1, nspins 446 rho_ao => rho_ao_kp(ispin, :) 447 CALL calculate_drho_elec(matrix_p_kp=rho_ao, & 448 drho=drho_r(3*(ispin - 1) + 1:3*ispin), & 449 drho_gspace=drho_g(3*(ispin - 1) + 1:3*ispin), & 450 qs_env=qs_env, soft_valid=gapw) 451 END DO 452 CALL qs_rho_set(rho_struct, drho_r_valid=.TRUE., drho_g_valid=.TRUE.) 453 ENDIF 454 455 ! if needed compute also the kinetic energy density 456 IF (dft_control%use_kinetic_energy_density) THEN 457 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp) 458 DO ispin = 1, nspins 459 rho_ao => rho_ao_kp(ispin, :) 460 CALL calculate_rho_elec(matrix_p_kp=rho_ao, & 461 rho=tau_r(ispin), & 462 rho_gspace=tau_g(ispin), & 463 total_rho=dum, & ! presumably not meaningful 464 ks_env=ks_env, soft_valid=gapw, & 465 compute_tau=.TRUE., & 466 task_list_external=task_list, & 467 pw_env_external=pw_env) 468 END DO 469 CALL qs_rho_set(rho_struct, tau_r_valid=.TRUE., tau_g_valid=.TRUE.) 470 ENDIF 471 END IF 472 473 ! GAPW o GAPW_XC require the calculation of hard and soft local densities 474 IF (gapw) THEN 475 CPASSERT(.NOT. PRESENT(task_list_external)) 476 CPASSERT(.NOT. PRESENT(pw_env_external)) 477 CALL qs_rho_get(rho_struct, rho_ao_kp=rho_ao_kp) 478 IF (PRESENT(local_rho_set)) THEN 479 CALL calculate_rho_atom_coeff(qs_env, rho_ao_kp, local_rho_set%rho_atom_set) 480 ELSE 481 CALL calculate_rho_atom_coeff(qs_env, rho_ao_kp) 482 ENDIF 483 ENDIF 484 IF (gapw_xc) THEN 485 CPASSERT(.NOT. PRESENT(task_list_external)) 486 CPASSERT(.NOT. PRESENT(pw_env_external)) 487 CALL get_qs_env(qs_env=qs_env, rho_xc=rho_xc) 488 CALL qs_rho_get(rho_xc, & 489 rho_ao_kp=rho_xc_ao, & 490 rho_r=rho_xc_r, & 491 rho_g=rho_xc_g, & 492 tot_rho_r=tot_rho_r_xc, & 493 drho_g=drho_xc_g, & 494 tau_r=tau_xc_r, & 495 tau_g=tau_xc_g) 496 CALL calculate_rho_atom_coeff(qs_env, rho_ao_kp) 497 ! copy rho_ao into rho_xc_ao 498 DO ispin = 1, nspins 499 DO img = 1, nimg 500 CALL dbcsr_copy(rho_xc_ao(ispin, img)%matrix, rho_ao_kp(ispin, img)%matrix) 501 END DO 502 END DO 503 DO ispin = 1, nspins 504 rho_ao => rho_xc_ao(ispin, :) 505 CALL calculate_rho_elec(matrix_p_kp=rho_ao, & 506 rho=rho_xc_r(ispin), & 507 rho_gspace=rho_xc_g(ispin), & 508 total_rho=tot_rho_r_xc(ispin), & 509 ks_env=ks_env, soft_valid=gapw_xc) 510 END DO 511 CALL qs_rho_set(rho_xc, rho_r_valid=.TRUE., rho_g_valid=.TRUE.) 512 ! if needed compute also the gradient of the density 513 IF (dft_control%drho_by_collocation) THEN 514 DO ispin = 1, nspins 515 rho_ao => rho_xc_ao(ispin, :) 516 CALL calculate_drho_elec(matrix_p_kp=rho_ao, & 517 drho=rho_xc_r(3*(ispin - 1) + 1:3*ispin), & 518 drho_gspace=drho_xc_g(3*(ispin - 1) + 1:3*ispin), & 519 qs_env=qs_env, soft_valid=gapw_xc) 520 END DO 521 CALL qs_rho_set(rho_xc, drho_r_valid=.TRUE., drho_g_valid=.TRUE.) 522 ENDIF 523 ! if needed compute also the kinetic energy density 524 IF (dft_control%use_kinetic_energy_density) THEN 525 DO ispin = 1, nspins 526 rho_ao => rho_xc_ao(ispin, :) 527 CALL calculate_rho_elec(matrix_p_kp=rho_ao, & 528 rho=tau_xc_r(ispin), & 529 rho_gspace=tau_xc_g(ispin), & 530 total_rho=dum, & ! presumably not meaningful 531 ks_env=ks_env, soft_valid=gapw_xc, & 532 compute_tau=.TRUE.) 533 END DO 534 CALL qs_rho_set(rho_xc, tau_r_valid=.TRUE., tau_g_valid=.TRUE.) 535 ENDIF 536 ENDIF 537 538 CALL timestop(handle) 539 540 END SUBROUTINE qs_rho_update_rho 541 542! ************************************************************************************************** 543!> \brief Duplicates a pointer physically 544!> \param rho_input The rho structure to be duplicated 545!> \param rho_output The duplicate rho structure 546!> \param qs_env The QS environment from which the auxiliary PW basis-set 547!> pool is taken 548!> \par History 549!> 07.2005 initial create [tdk] 550!> \author Thomas D. Kuehne (tkuehne@phys.chem.ethz.ch) 551!> \note 552!> Associated pointers are deallocated, nullified pointers are NOT accepted! 553! ************************************************************************************************** 554 SUBROUTINE duplicate_rho_type(rho_input, rho_output, qs_env) 555 556 TYPE(qs_rho_type), POINTER :: rho_input, rho_output 557 TYPE(qs_environment_type), POINTER :: qs_env 558 559 CHARACTER(len=*), PARAMETER :: routineN = 'duplicate_rho_type', & 560 routineP = moduleN//':'//routineN 561 562 INTEGER :: handle, i, nspins, rebuild_each_in 563 LOGICAL :: drho_g_valid_in, drho_r_valid_in, rho_g_valid_in, rho_r_valid_in, soft_valid_in, & 564 tau_g_valid_in, tau_r_valid_in 565 REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_g_in, tot_rho_g_out, & 566 tot_rho_r_in, tot_rho_r_out 567 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_in, rho_ao_out 568 TYPE(dft_control_type), POINTER :: dft_control 569 TYPE(pw_env_type), POINTER :: pw_env 570 TYPE(pw_p_type), DIMENSION(:), POINTER :: drho_g_in, drho_g_out, drho_r_in, drho_r_out, & 571 rho_g_in, rho_g_out, rho_r_in, rho_r_out, tau_g_in, tau_g_out, tau_r_in, tau_r_out 572 TYPE(pw_p_type), POINTER :: rho_r_sccs_in, rho_r_sccs_out 573 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 574 575 CALL timeset(routineN, handle) 576 577 NULLIFY (dft_control, pw_env, auxbas_pw_pool) 578 NULLIFY (rho_ao_in, rho_ao_out) 579 NULLIFY (rho_r_in, rho_r_out, rho_g_in, rho_g_out, drho_r_in, drho_r_out) 580 NULLIFY (drho_g_in, drho_g_out, tau_r_in, tau_r_out, tau_g_in, tau_g_out) 581 NULLIFY (tot_rho_r_in, tot_rho_r_out, tot_rho_g_in, tot_rho_g_out) 582 NULLIFY (rho_r_sccs_in, rho_r_sccs_out) 583 584 CPASSERT(ASSOCIATED(rho_input)) 585 CPASSERT(ASSOCIATED(rho_output)) 586 CPASSERT(ASSOCIATED(qs_env)) 587 CPASSERT(qs_env%ref_count > 0) 588 589 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, dft_control=dft_control) 590 CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool) 591 nspins = dft_control%nspins 592 593 CALL qs_rho_clear(rho_output) 594 595 CALL qs_rho_get(rho_input, & 596 rho_ao=rho_ao_in, & 597 rho_r=rho_r_in, & 598 rho_g=rho_g_in, & 599 drho_r=drho_r_in, & 600 drho_g=drho_g_in, & 601 tau_r=tau_r_in, & 602 tau_g=tau_g_in, & 603 tot_rho_r=tot_rho_r_in, & 604 tot_rho_g=tot_rho_g_in, & 605 rho_g_valid=rho_g_valid_in, & 606 rho_r_valid=rho_r_valid_in, & 607 drho_g_valid=drho_g_valid_in, & 608 drho_r_valid=drho_r_valid_in, & 609 tau_r_valid=tau_r_valid_in, & 610 tau_g_valid=tau_g_valid_in, & 611 rho_r_sccs=rho_r_sccs_in, & 612 soft_valid=soft_valid_in, & 613 rebuild_each=rebuild_each_in) 614 615 ! rho_ao 616 IF (ASSOCIATED(rho_ao_in)) THEN 617 CALL dbcsr_allocate_matrix_set(rho_ao_out, nspins) 618 CALL qs_rho_set(rho_output, rho_ao=rho_ao_out) 619 DO i = 1, nspins 620 ALLOCATE (rho_ao_out(i)%matrix) 621 CALL dbcsr_copy(rho_ao_out(i)%matrix, rho_ao_in(i)%matrix, & 622 name="myDensityMatrix_for_Spin_"//TRIM(ADJUSTL(cp_to_string(i)))) 623 CALL dbcsr_set(rho_ao_out(i)%matrix, 0.0_dp) 624 END DO 625 END IF 626 627 ! rho_r 628 IF (ASSOCIATED(rho_r_in)) THEN 629 ALLOCATE (rho_r_out(nspins)) 630 CALL qs_rho_set(rho_output, rho_r=rho_r_out) 631 DO i = 1, nspins 632 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r_out(i)%pw, & 633 use_data=REALDATA3D, in_space=REALSPACE) 634 rho_r_out(i)%pw%cr3d(:, :, :) = rho_r_in(i)%pw%cr3d(:, :, :) 635 END DO 636 END IF 637 638 ! rho_g 639 IF (ASSOCIATED(rho_g_in)) THEN 640 ALLOCATE (rho_g_out(nspins)) 641 CALL qs_rho_set(rho_output, rho_g=rho_g_out) 642 DO i = 1, nspins 643 CALL pw_pool_create_pw(auxbas_pw_pool, rho_g_out(i)%pw, & 644 use_data=COMPLEXDATA1D, & 645 in_space=RECIPROCALSPACE) 646 rho_g_out(i)%pw%cc(:) = rho_g_in(i)%pw%cc(:) 647 END DO 648 END IF 649 650 ! SCCS 651 IF (ASSOCIATED(rho_r_sccs_in)) THEN 652 CALL qs_rho_set(rho_output, rho_r_sccs=rho_r_sccs_out) 653 CALL pw_pool_create_pw(auxbas_pw_pool, rho_r_sccs_out%pw, & 654 in_space=REALSPACE, & 655 use_data=REALDATA3D) 656 rho_r_sccs_out%pw%cr3d(:, :, :) = rho_r_sccs_in%pw%cr3d(:, :, :) 657 END IF 658 659 ! drho_r and drho_g are only needed if calculated by collocation 660 IF (dft_control%drho_by_collocation) THEN 661 ! drho_r 662 IF (ASSOCIATED(drho_r_in)) THEN 663 ALLOCATE (drho_r_out(3*nspins)) 664 CALL qs_rho_set(rho_output, drho_r=drho_r_out) 665 DO i = 1, 3*nspins 666 CALL pw_pool_create_pw(auxbas_pw_pool, drho_r_out(i)%pw, & 667 use_data=REALDATA3D, in_space=REALSPACE) 668 drho_r_out(i)%pw%cr3d(:, :, :) = drho_r_in(i)%pw%cr3d(:, :, :) 669 END DO 670 END IF 671 672 ! drho_g 673 IF (ASSOCIATED(drho_g_in)) THEN 674 ALLOCATE (drho_g_out(3*nspins)) 675 CALL qs_rho_set(rho_output, drho_g=drho_g_out) 676 DO i = 1, 3*nspins 677 CALL pw_pool_create_pw(auxbas_pw_pool, drho_g_out(i)%pw, & 678 use_data=COMPLEXDATA1D, & 679 in_space=RECIPROCALSPACE) 680 drho_g_out(i)%pw%cc(:) = drho_g_in(i)%pw%cc(:) 681 END DO 682 END IF 683 END IF 684 685 ! tau_r and tau_g are only needed in the case of Meta-GGA XC-functionals 686 ! are used. Therefore they are only allocated if 687 ! dft_control%use_kinetic_energy_density is true 688 IF (dft_control%use_kinetic_energy_density) THEN 689 ! tau_r 690 IF (ASSOCIATED(tau_r_in)) THEN 691 ALLOCATE (tau_r_out(nspins)) 692 CALL qs_rho_set(rho_output, tau_r=tau_r_out) 693 DO i = 1, nspins 694 CALL pw_pool_create_pw(auxbas_pw_pool, tau_r_out(i)%pw, & 695 use_data=REALDATA3D, in_space=REALSPACE) 696 tau_r_out(i)%pw%cr3d(:, :, :) = tau_r_in(i)%pw%cr3d(:, :, :) 697 END DO 698 END IF 699 700 ! tau_g 701 IF (ASSOCIATED(tau_g_in)) THEN 702 ALLOCATE (tau_g_out(nspins)) 703 CALL qs_rho_set(rho_output, tau_g=tau_g_out) 704 DO i = 1, nspins 705 CALL pw_pool_create_pw(auxbas_pw_pool, tau_g_out(i)%pw, & 706 use_data=COMPLEXDATA1D, & 707 in_space=RECIPROCALSPACE) 708 tau_g_out(i)%pw%cc(:) = tau_g_in(i)%pw%cc(:) 709 END DO 710 END IF 711 END IF 712 713 CALL qs_rho_set(rho_output, & 714 rho_g_valid=rho_g_valid_in, & 715 rho_r_valid=rho_r_valid_in, & 716 drho_g_valid=drho_g_valid_in, & 717 drho_r_valid=drho_r_valid_in, & 718 tau_r_valid=tau_r_valid_in, & 719 tau_g_valid=tau_g_valid_in, & 720 soft_valid=soft_valid_in, & 721 rebuild_each=rebuild_each_in) 722 723 ! tot_rho_r 724 IF (ASSOCIATED(tot_rho_r_in)) THEN 725 ALLOCATE (tot_rho_r_out(nspins)) 726 CALL qs_rho_set(rho_output, tot_rho_r=tot_rho_r_out) 727 DO i = 1, nspins 728 tot_rho_r_out(i) = tot_rho_r_in(i) 729 END DO 730 END IF 731 732 ! tot_rho_g 733 IF (ASSOCIATED(tot_rho_g_in)) THEN 734 ALLOCATE (tot_rho_g_out(nspins)) 735 CALL qs_rho_set(rho_output, tot_rho_g=tot_rho_g_out) 736 DO i = 1, nspins 737 tot_rho_g_out(i) = tot_rho_g_in(i) 738 END DO 739 740 END IF 741 742 CALL timestop(handle) 743 744 END SUBROUTINE duplicate_rho_type 745 746END MODULE qs_rho_methods 747